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Abstract 


Implementation  of  three  algorithms,  (i)  moment  equation,  (ii)  Monte  Carlo,  and  (iii) 
quantum  Liouville  equation  algorithms,  were  used  in  a  program  to  determine  the  high 
speed  and  high  frequency  operation  of  submicron  electron  devices.  For  a  pseudomorphic 
high  electron  mobility  transistor,  high  frequency,  small  signal,  subpicosecond  charge  density 
waves  were  observed  to  form  within  the  two  dimensional  electron  gas.  Large  signal 
operation  of  the  PHEMT  indicated  that  the  switching  time  of  the  device  was  governed  by 
the  longest  relaxation  effect,  the  energy  relaxation  time,  estimated  to  be  longer  than  two 
picoseconds.  A  simple  two  terminal  device  configuration  was  examined.  It  was  determined 
that  measurements  of  its  transient  behavior,  would  expose  differences  in  the  key  relaxation 
times  governing  III-V  device  behavior,  and  provide  the  first  direct  measurement  of 
nonequilibrium  effects  in  semiconductor  devices. 
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Numerical  Simulation  of  the  Function  of  Scientific 

« 

Instrumentation  for  Measuring  the  Speed  of  Electron  Devices 

1 .  Introduction 

The  development  of  techniques  to  generate  laser  pulses  with  picosecond  and 
sub-picosecond  duration  has  finally  provided  the  means  to  probe  ultrafast  relaxation 
process  in  semiconductors  to  determine  the  response  time  of  new  high  speed  devices. 
This  process  has  introduced  a  series  of  significant  dilemmas.  It  was  thought  that  as  device 
size  decreases  the  response  time  decreases  and  the  speed  of  the  device  increases.  Effects 
such  as  velocity  overshoot  entered  prominently.  But  we  are  instead  finding  several  new 
issues.  It  is  possible  to  reduce  the  response  time  of  devices;  but  we  are  finding  that  when 
very  small  size  devices  enter  the  picture,  capacitive  and  tunnelling  times  conspire  to  make 
tiny  feature  size  devices  slow;  that  is  response  times  are  often  greater  than  one 
picosecond,  although  higher  frequency  components  do  exist. 

Thus  the  question  that  remains  for  all  devices  is:  how  fast  can  we  transport  charge  from 
one  end  of  the  structure  to  the  next,  and  how  much  charge  can  be  transported ?  In  an  earlier 
study  (Phase  I  SBIR)  it  was  suggested  that  the  shortest  response  time  of  a  device  was  the 
time  it  takes  a  device  to  approach  5-10  %  of  its  steady  state  value.  This  is  the  response 
time  for  large  signal  behavior,  but  it  is  not  necessarily  the  response  time  for  small  signal 
behavior.  Indeed  carriers  may  be  expected  to  follow  almost  any  signal;  but  whether  there 
will  be  any  gain  or  any  useful  operation  is  to  be  determined. 

To  determine  the  speed  signature  of  a  device  we  must  understand  those  features  of  device 
operation  that  affect  its  transient  behavior.  This  was  the  broad  motivation  for  the  Phase  II 
SBIR  study  that  this  report  summarizes.  In  particular  through  three  distinctly  different 
algorithms:  (a)  moment  algorithms,  (b)  Monte  Carlo  algorithms,  and  (c)  quantum 
transport  algorithms,  we  examined  the  transient  device  behavior  of  the  pseudomorphic 
HEMT,  and  were  able  to  suggest  a  direct  experimental  means  of  identifying 
nonequilibrium  transient  effects  in  semiconductors.  It  is  worthwhile  dwelling  on  this  last 
point. 

For  many  years  the  device  and  scientific  community  has  been  seeking  a  measure  of 
velocity  overshoot.  The  conundrum  is  that  most  think  that  velocity  overshoot  exists,  and 
that  the  classical  steady  state  description  of  device  behavior  underestimates  velocity  and 
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hence  current.  But  there  was  no  way  to  compare  the  two  because  only  one  velocity  was 
correct.  Thus  if  velocity  overshoot  always  exists  how  do  we  detect  it ?  Clearly,  the  speed  of 
the  device  is  not  a  measure  of  velocity  overshoot,  since  in  general  current  is  measured,  not 
velocity.  To  address  this  we  have  taken  a  different  tack  and  suggest,  as  a  result  of  this 
study,  that,  rather  than  attempt  to  measure  velocity  overshoot,  measurements  should 
concentrate  on  determining  the  differences  between  energy  and  momentum  relaxation 
times,  since  these  differences,  when  significant,  are  responsible  for  nonequilibrium  device 
behavior. 

But  in  addition  to  providing  a  means  by  which  nonequilibrium  effects  i  ccur,  significant 
code  development  pointed  to  the  existence  of  high  frequency  charge  density  waves  in  the 
channel  of  the  pseudomorphic  HEMT.  These  waves  corresponded  to  local  density 
changes  within  the  two  dimensional  electron  gas,  and  suggest  high  frequency  terahertz 
response.  The  behavior  of  these  oscillation  also  suggest  strong  high  frequency  capacitive 
contributions,  but  the  details  indicate  otherwise. 

In  the  course  of  the  study  both  Monte  Carlo  and  moment  equation  algorithms  were 
examined.  Steady  state  curves  of  the  two  were  obtained  for  simple  (non  heterostructure 
FETs)  and  it  was  found  that  virtually  identical  results  could  be  obtained.  Because 
transient  nonequilibrium  codes  were  robust,  the  transient  operation  was  examined  using 
the  moment  codes.  The  primary  three  terminal  structure  studied  during  the  Phase  II 
program  involved  a  two  dimensional  electron  gas,  and  another  issue  was  how  do  quantum 
contributions  affect  the  distribution  of  charge ?.  Originally,  this  problem  was  to  be 
addressed  through  incorporating  quantum  effects  into  the  scattering  integrals.  While  this 
was  done,  most  of  the  effort  was  devoted  to  solving  the  multiparticle  quantum  Liouville 
equation,  which  is  a  generalization  of  the  Schrodinger’s  equation,  for  the  heterostructure 
configuration  associated  with  the  pseudo-morphic  HEMT. 

One  of  the  key  issues  addressed  during  the  study  was  whether  nonequilibrium  transients 
could  be  measured.  As  a  result  of  the  simulations  it  was  determined  that  it  is  possible  to 
expose  the  different  relaxation  times  in  a  common  device  configuration  in  a  manner  that 
would  lead  to  the  first  experimental  demonstration  of  nonequilibrium  contributions  to 
velocity  overshoot. 

The  basic  conclusion  of  the  Phase  I  study  is  reaffirmed  for  the  Phase  II  study;  namely  a 
measure  of  a  devices  speed  is  the  time  it  takes  a  device  to  approach  5-10  %  of  its  steady 
state  value. 
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2.  The  Governing  Equations 

The  Moment  Equations 


Introduction 

As  indicated  in  our  introductory  remarks  the  thrust  of  this  study  was  in  both  the  steady 
state  and  transient  behavior  of  devices.  The  equations  governing  this  transport  must  be 
relevant.  In  choosing  the  equations  care  was  taken  with  respect  to  the  device  physics,  and 
computational  resources.  The  three  terminal  structure  that  we  are  dealing  with  is  a 
pseudomorphic  HEMT  whose  operation  requires  quantum  transport  for  it’s  description. 
In  earlier  studies,  moments  of  the  Wigner  function  were  introduced  and  the  effects  of 
quantum  transport  were  introduced  through  the  quantum  potential  (Grubin  and 
Kreskovsky,  1989).  The  use  of  the  quantum  potential  has  not  progressed  far  enough  and 
was  not  implemented  during  the  study,  although  quantum  transport  calculations  were 
performed  through  the  incorporation  of  the  density  matrix.  These  quantum  transport 
studies  were  introduced  primarily  to  ascertain  the  significance  of  some  of  the  space 
charge  profiles  obtained  through  the  more  traditional  means  of  transport,  and  regarded  as 
supplementary  to  the  core  simulations. 

The  primary  transport  calculations  as  performed  today  are  classical  calculations  in  which 
quantum  transport  is  incorporated  into  the  scattering  integrals,  and  into  the  treatment  of 
carriers  in  different  subbands  as  interacting  though  scattering  process.  In  such  situations  a 
carrier  traversing  a  region  of  a  device  in  which  it  makes  a  transition  from  a  wide  band  gap 
region  to  a  narrow  band  gap  region  is  thought  to  have  all  of  its  bandgap  energy  converted 
into  kinetic  energy.  Detailed  quantum  transport  calculations  that  incorporate  the 
quantum  potential  indicate  that  this  is  not  occurring.  Instead  there  are  quantum 
mechanical  forces  that  exist  on  either  side  of  the  wide  band  gap  region  that  tend  to 
mitigate  the  transition  in  energy  across  the  interface.  These  effects  are  not  included  in  the 
Monte  Carlo  calculations  nor  are  they  included  in  those  moment  equations  that  do  not 
include  quantum  corrections.  But  for  the  most  part  they  are  all  that  we  have  to  deal  with 
at  the  present  time.  And  so  the  quantum  calculations  are  introduced  to  determine  how 
well  the  classical  transport  calculations  represent  some  of  the  quantum  streaming 
contributions 

During  the  course  of  this  study  the  Monte  Carlo  algorithm  was  implemented  and  applied 
to  a  simple  FET,  as  discussed  in  the  introduction.  This  same  FET  was  examined  using  the 
moments  of  the  Boltzmann  transport  equation.  The  results  of  the  simulation  indicated 
that  both  results  would  yield  substantially  similar  results  quantitatively;  the  results  were 
certainly  qualitatively  similar.  Insofar  as  much  of  the  one  dimensional  calculations 
indicate  better  numerical  control  in  representing  boundary  conditions,  it  was  decided  to 
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devote  most  of  the  simulation  effort  to  enhance  the  moment  algorithms,  in  particular  with 
regard  to  analyzing  short  time  phenomena  on  a  scale  not  represented  in  the  current 
literature.  It  was  also  felt  that  this  exercise  would  bear  further  advantages  in  that  the 
moment  algorithms  are  currently  capable  of  modification  to  include  quantum  corrections, 
whereas  the  Monte  Carlo  algorithms  are  not.  Thus  further  development  of  the  moment 
algorithms  would  position  us  to  address  more  realistic  problems  on  a  shorter  time  scale. 

The  moment  algorithm  used  during  the  study  was  the  three  dimensional  moment 
algorithm  developed  at  SRA,  but  modified  to  deal  with  the  transients  in  three  terminal 
devices  in  a  more  accurate  manner.  During  the  course  of  the  study  this  algorithm  was 
generalized  to  deal  with  two  dimensional  transport  as  well  as  three  dimensional  transport. 
The  two  dimensional  transport  algorithm  was  not  implemented  during  the  course  of  the 
study,  but  it  is  anticipated  that  it  will  be  used  heavily  in  conjunction  with  quantum 
corrective  transport.  The  development  of  the  two  dimensional  transport  as  well  as  the 
three  dimensional  equations  are  discussed  in  this  section.  Furthermore,  while  we  have 
developed  the  equations  for  nonparabolic,  nonspherical  systems,  the  discussion  below  is 
confined  to  parabolic  systems.  The  new  work  discussed  below  is  for  two  dimensional 
transport.  Here  there  are  two  elements  to  consider;  first,  the  streaming  terms,  and  second , 
the  scattering  integrals.  The  streaming  terms  have  been  discussed  many  times  by  many 
authors,  and  are  included  here  for  completeness.  The  scattering  integrals  are  new.  By 
new  we  mean  the  following.  There  has  been  a  considerable  effort  undertaken  to  include 
scattering  in  the  two  dimensional  electron  gas  for  incorporation  into  Monte  Carlo 
algorithms,  and  some  of  that  work  was  undertaken  during  this  study.  This  is  a 
straightforward  prescription  and  simply  serves  to  enhance  a  numerical  capability.  Indeed 
in  many  hydrodynamic  transport  models  in  semiconductors  Monte  Carlo  algorithms  have 
been  used  to  obtain  the  scattering  integrals  relevant  to  transport.  At  SRA  a  different 
approach  has  always  been  taken.  From  the  very  beginning  of  the  SRA  studies  with  the 
moments  of  the  Boltzmann  transport  equation,  scattering  integrals  were  computed 
assuming  a  displaced  Maxwellian.  In  every  case  when  a  comparison  was  made  to  the 
Monte  Carlo  calculations,  minor  adjustments  in  the  input  parameters  such  as  the 
intervalley  coupling  coefficient  was  all  that  was  needed  to  achieve  similar  results  from  the 
two  sets  of  calculations.  It  was  with  this  confidence  that  a  set  of  scattering  integrals  was 
developed.  The  analysis  discussed  is  for  one  type  of  scattering  event,  the  polar  phonon. 

The  discussion  below  is  geared  to  two  dimensional  transport;  but  for  assurances  that  the 
procedures  were  correct,  parallel  three  dimensional  transport  arguments  were  used  as 
check  points.  In  the  discussion  both  two  and  three  dimensional  transport  terms  are 
discussed.  Note  all  three  dimensional  contributions  are  part  of  the  current  SRA  moment 
equation  algorithm. 
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The  Boltzmann  Transport  Equation  and  it’s  Moments 

All  information  is  obtained  from  the  Boltzmann  tranport  equation: 

(2.1)  .  {a/at  +  [ftk/m*]- divx +[qF/ft]*  divjc}f(x,k,t)  =  a[f(x,k,t)]/at|cop 

The  moments  of  the  BTE  are  obtained  by  multiplying  equation  (2.1)  by  successive  powers 
of  k,  and  integrating  over  k-space. 

(2.2)  a  {nj  <  <t>(* )  >  i}/8 1 +vr{nj  <  (v^E^C*  )fti  >  j} 

=  qFM-nj<Vj^(A)>j+  <anj3  <+>i(l)/dt>00n 
For  three  dimensional  transport: 

(2.3)  nj  <  *(* )  >  i  =  [2/(2*) 3  ]J<p(A  )(k)fi(k)d3ki 
and  for  two  dimensional  systems: 

(2.4)  ni<*(*)>i  =  [2/(2it)2  U  <t>(*  )(k)fj[(k)d2  kj 
For  parabolic  bands: 

(3D)  d3k  =  sinfld0d^k2dk 

(2D)  d2  k  =  da  kdk 

In  three  dimensions  the  density  is  per  cubic  centimeter;  in  two  dimension  the  density  is  per 
square  centimeter.  As  an  example  consider  an  evaluation  of  the  Fermi  momentum  at 
absolute  zero  for  two  and  three  dimensional  systems  (here  we  deal  with  <t>(° )  =  1).  For 
both  two  and  three  dimensional  systems  the  distribution  is  equal  to  unity  up  to  energies 
equal  to  the  Fermi  energy,  and  zero  beyond.  For  parabolic  bands  the  integration  is  a 
volume  integration  with  a  value  equal  to  [4/3x]kf3 ,  whereas  for  the  circular  symmetry  in  the 
case  of  the  two  dimensional  systems  the  integral  is  equal  to  xkf2.  Thus  for  three 
dimensional  systems: 

(2.5)  kf  =  [3ir2ni(3)]1/3 
and  for  two  dimensional  systems: 

(2.6)  kf  =  [2xni(2)]1/2 
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The  scattering  integrals  for  three  dimensional  systems  are 

(2.7)  <  d  nj  <  <t>  >  {(*■  )/d  t  >  con  = 

[l/4r>  ];<,(<  )d=  kj[V/8w3]/d3  kj{  W(ki,kj>f(kj)[l-f(ki)]  -  W(kj,kj)f(ki)[l-f(kj))} 
and  for  two  dimensional  systems: 

(2.8)  <3n|<0>j(^)/3t>coi|  = 

[1/2*  » !/*(■*  )d’  kjt/l/4*  » ]/d5kj{W(kj,kj)f(kj)[l-f(kj)]  -  W(kj,ki)f(ki)[l-f(kj)]} 

In  the  above  the  indices  i  and  j  account  for  intervalley  scattering.  For  i=j,  there  is 
intravalley  scattering. 

For  the  purpose  of  vector  direction,  we  will  assume  that  a  two  dimensional  electron  gas  has 
formed.  That  the  coordinate  V  represents  transport  normal  to  the  gate,  that  V  represents 
transport  in  the  source  to  the  drain  contact,  and  that  ’y’  represents  transport  in  the 
unmodelled  direction.  The  Fermi  energy  for  the  two  dimensional  system  defines  the  number 
of  states  per  square  centimeter . 

In  the  absence  of  collisions,  an  exact  solution  to  the  BTE  can  be  obtained  by  solving  an 
infinite  set  of  coupled  moment  equations.  This  procedure  is  now  standard  and  has  been 
discussed  in  a  variety  of  text  (see,  e.g.  Chapman  and  Cowling,  1964).  This  method  is  not, 
however,  satisfactory  and  a  variety  of  sophisticated  techniques  are  being  developed, 
particularly  Monte  Carlo  procedures  (see,  e.g.,  Jacoboni  and  Reggiani,  1985),  are  being 
invoked.  The  procedure  used  here  is  based  on  a  philosophy  discussed  in  earlier  papers  by 
the  author  (see,  e.g.  Grubin  et.  al.  (1982)  and  (1985))  and  involved  first  estimating  the 
distribution  function  and  then  taking  moments.  The  specific  form  of  the  moment 
equations  are  sensitive  to  the  approximations  used  for  estimating  the  distribution  function, 
as  will  be  apparent  below.  The  distribution  function  most  often  used  in  examining 
transient  nonlinear  high  field  transport  is  the  displaced  Maxwellian.  This  is  the  approach 
taken  below. 

The  Displaced  Maxwellian  and  the  Moment  Equations 

In  calculating  the  transient  properties  of  submicron  devices  we  assumed  that  the 
nonequilibrium  distribum  function  is  of  the  form  of  a  displaced  Maxwellian: 

(2.9)  folE(k-kd)]  =  (exp-Py  -  gradk])f0(E) 
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where 


(2.10)  f0(E)  =  Aexp-(E/kBTe) 

m 

where  A  is  a  constant  to  be  determined  and  is  different  for  both  two  and  three  dimensional 
systems.  Te  is  an  electron  temperature.  For  parabolic  spherical  energy  bands  in  three 
dimensions: 

(2.11a)  E(k)  =  [«2/2m](kx2+ky2+kz2) 

In  two  dimensions: 

(2.11b)  E(k)  =  [ft  2/2m](kx2  +  kz2 ) 

The  first  integral  of  interest  involves  the  determination  of  the  contant  A.  Its  evaluation, 
which  is  equivalent  to  determining  the  form  of  the  first  moment  (see  Eq.  (4))  is  dependent 
on  the  requirement  that 

(2.12)  -  N 

where  N  is  the  total  number  of  particles  under  study.  Note,  we  have  suppressed  the  energy 
notation  in  the  distribution  function. 

For  parabolic  surfaces 

(2.13)  N  =Xkf„(k-kd)  =  [2VA(3  )/( 2„)3  ]Jexp-[E(k-k<|)/kBTe]d3  k 

(2.14)  N  =Xkf0(k-kd)  =  [2AA(*  )/<2»)>  ]Jexp-[E(k-kd)/kBTe]d  =  k 

where  the  factor  "2"  arises  because  each  quantum  state  "k"  can  hold  two  electrons  of 
different  spin.  We  now  defined  three  and  two  dimensional  densities: 

(2. 15a)  n(3)  =  N/V  =  [A(s  )/4* 3  ]/d 3  k  exp-[E(k-kd)/kBTe] 

where  n  =  N/V  is  the  density  of  carriers,  and  for  two  dimensional  systems: 

(2. 16b)  n(2)  =  N/A  =  [a(  2  )/2it 2  ]Jd 2  k  exp-[E(k-kd)/kBTe] 

where  for  a  three  dimensional  system: 


-7- 


(2.17) 


d3k  =  sin6d6d'ir[2m/ft2]3/2  JEdE/2 


For  a  two  dimensional  system 
(2.18)  d2k  =  de[m/ft  2]dE 

Thus  with  the  introduction  of  the  thermal  deBroglie  wavelength: 


(2.19) 

X2  =  ftVCmdkBTe) 

(2.20) 

A(3)  =  [4W2][X3]n(3) 

and 

(2.21) 

A(2)  =  [2*][X2]n(2) 

Having  established  the  value  of  the  coefficients  for  the  two  and  three  dimensional 
systems  it  is  direct  to  demonstrate  that  for  parabolic  systems  the  mean  current  is  given 
by  neft  k^/m,  and  that  the  pressure  gradient  for  two  and  three  dimensional  systems  is 
the  same.  Thus  the  first  two  moment  equations  are  identical  in  form  with  the 
exception  of  the  collision  integrals. 

Thus  the  first  moment  equation,  for  carrier  balance,  is: 

(2.22)  a  nj/a  t-divj  j/e  =  [a  n  [la  t]cou 

The  second  moment  equation,  for  momentum  balance,  is: 

(2.23)  [a/at  +  (div«vj)]nj«k(ji  =  qnjF  -  div  njkgTj  +  [anjftkj/at]cou 
For  energy  balance,  and  for  a  displaced  Maxwellian,  the  three  dimensional  energy  is: 

(2.24)  Wj  =  [3/2]njkBTe  +  n,ft 2 k^7  /2m 
For  two  dimensional  transport: 

(2.25)  Wj  =  njkBTe  +  ni«2ktj2/2m 

A  result  that  is  expected  based  upon  equipartition.  Thus  the  energy  balance  equations  are: 
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(2.26) 


[d/3t  +  (div*vj)]Wj  =  qnjvj*  F  -vj*  divnjkgTj  +  [3  Wj/3t]co]] 


with  the  appropriate  quantity  identifying  energy. 

The  key  is^ue  in  evaluating  the  moment  equations  are  the  scattering  integrals,  and  this 
represents  a  new  effort  not  discussed  in  the  literature,  as  most  discussion  are  confined  to 
using  Monte  Carlo  evaluations  of  the  scattering  integrals  in  the  moment  equations.  We 
have  taken  the  Butcher  (1967)  approach  and  have  evaluated  the  three  dimensional 
integrals  and  have  begun  to  evaluate  the  two  dimensional  integu  The  evaluation  is 
tedious  and  only  a  brief  disucssion  is  given  here. 

The  Scattering  Integrals 

Here  we  will  concentrate  only  on  the  momentum  balance  scattering  integral.  We  define  the 
momentum  scattering  rate  implicitly  as : 

(2.27)  k»X  ft (k’-k)  W(k’,k)  =  -ftkr(E) 

where  r(E)  denotes  a  scattering  rate  and  is  energy  dependent.  The  specific  form  of  the 
energy  dependence  is  dependent  upon  the  scattering  mechanism.  The  quantity  r  (E)  is  not 
generally  the  inverse  of  the  mean  free  time  between  collisions.  It  is  the  rate  at  which  the  net 
momentum  is  destroyed  as  a  result  of  scattering.  The  scattering  rate  is  used  in  evaluating 
the  field  dependent  mobility. 

With  the  scattering  rate  defined  as  in  Eq.  (2.27),  the  mean  time  rate  of  change  of 
momentum  due  to  collisions  is  given  by 

(2.28)  [d  njf.  k  j/a  t]coll  =  -[1/(4* 3  )]Jd3  kft  k T  (E)f(k) 

For  two  dimensional  transport: 

(2.29)  (a  njdkj/a  Dell = -[1/(2* *  )lfd*  k*  kr(E)f(k) 

Note  for  r(E)  equal  to  a  constant  1/r,  independent  of  energy: 

(2.30)  [d  njft  kj/3  t]con  =  -njfi  ly/r 
For  three  dimensions,  defining: 

(2.31)  U(E)  =  3(f  coshr  -sinhC  )/f 3 
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where 


(2.32)  C  =*ld[(2raiE)*/*)/[mikBTe] 
and  with  x  =  E/kgTe: 

(2.33)  f=2Xkd(Jx) 
the  scattering  takes  the  form 
(2.34a) 

[a  njfi  kj/a  t]coxi  =  -njft  1/T  (5/2)]  [exp-[X 2  k^  2  ]  0  /“  dxT  (x)x3  /2  U(x)exp-x 

This  is  the  generic  form  of  the  scattering  integral  used  in  the  SRA  algorithms  for  the  past 
decade.  In  the  limiting  situation  as  kd=>0,  UW1  and  equation  (2.34a)  simplifies  to 

(2.34b)  [a  njft  kj/a  t]con  =  -njft  k^[  1/T  (5/2)]  0  J®  dxT  (x)x3  /2  exp-x 

For  two  dimensions  with  the  following  definition: 

(2.35)  Ij(f )  =  [l/2ic]0  J  2  *dz[exp + (f  cosz)]cosz 

where  Ij(f  )  is  a  modified  Bessel  function  of  order  T  (Abramowitz  and  Stegun,  1964)  we 
rewrite  the  scattering  term  as: 

(2.36) 


[a  nj«  kj/a  t]con = -[2nX  2  ]exp-[X 2  ly 2  ]i0 /kdk(*  kr  (E))Ix(f  )exp-[E/kBTe] 

Equation  (2.36)  is  a  new  expression  for  scattering.  We  rewrite  equation  (2.36  in  a  different 
form,  recognizing  that  in  the  limit  as  f=»0,  Ii(f  )=*f/2.  We  define  X(f  )-2Ij(f  )/f .  Then: 

(2.37a)  [anj*kj/at]cou  =  -i0n*kdexp-[X  2kd2]0/®dxxT(x)X(x)exp-x 

To  find  an  expression  for  the  low  field  mobility  we  are  interesting  in  evaluating  the  above 
in  the  limit  as  kd»0,  or  f  =✓  0,  X(f  )  =>  1.  Then: 

(2.37b)  [a  njft  kj/a  t]cou  =  -i0rtf»  k^0  J®  dxxT  (x)exp-x 

Thus  to  evaluate  the  above  integrals  for  use  in  the  moment  equations  we  need  an 
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expression  for  r(E).  For  the  case  of  three  dimensional  transport  and  polar  phonons  r(E)  is 
in  textbooks  (see,  e.g.,  Ridley,  1982)  and  we  simply  write  down  the  result: 


(2.38)  r(0(E)  =  Fo/[2mE]*[Nq{*  +  -(*«</E)sinh'1  (EMw0)*} 

+  (Nq+l){*-+  (ft£J0/E)sinh_1  (E/*w0-l)^  }] 

where: 

(2.39)  X±  =  (1±IM0I /E/- 
and  the  Frohlich  force  is  defined  as: 

(2.40)  F0  =  [ftwo/aoK,c«'1’/co‘1) 

The  scattering  of  two  dimensional  electrons  is  considered  next. 

Scattering  of  two  dimensional  electrons  by  three  dimensional  phonons 

As  in  the  case  of  scattering  of  three  dimensional  electrons  by  three  dimensional  phonons, 
the  situation  of  scattering  of  two  dimensional  electrons  by  three  dimensional  phonons 
proceeds  along  the  same  path.  As  discussed  by  Yokoyama  and  Hess  (1986),  and  by  Das 
Sarma  et  al  (1988),  the  matrix  element  for  scattering  two  dimensional  electrons  by  three 
dimensional  electrons  involves  a  sum  over  all  phonon  components  recognizing  that  the 
selection  rules  are  only  for  those  components  of  the  phonon  wave  vector  that  are  in  the 
plane  of  the  two  dimensional  electrons.  How  do  we  treat  this  and  effectively  define  a 
transition  probability  that  represents  two  dimensional  transport?  We  recognize  that  while 
we  are  interested  in  in-plane  transport,  if  we  perform  a  sum  over  all  phonon  wave  vectors 
then  the  phonon  sum  is  a  three  dimensional  sum.  Thus  we  are  interested  in  evaluating: 

(2.41)  k£*  (k’-k)W(k’,k)  =  [V/8*  3  )ffl  (k-k’)W(k,k’)dqyd2  Q 
where  Q  denotes  the  in-plane  phonon  wave  vector  with  components: 

Q  =  (Qx’Qz) 

Suitable  selection  rules  will  be  discussed  below  for  the  in-plane  k  vectors  and  the  in-plane 
phonon  vectors.  We  then  rewrite  equation  (2.41)  as 

(2.42)  k£  ft  (k’-k)W(k’,k)  =  [. AJ4 *  *  ]J«  (k’-k)w(k,k’)d 2  Q 
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where 


(2.43)  w(k,k’)  =  [l/2n][VM]/  dqyW(k,k’) 

If  we  include  a  suitable  overlap  integral  |  Imn(tlz)  1 2 » the  matrix  element  for  scattering  the 
two  dimensional  electrons  is  effectively: 

(2.44)  w(k,k’)  =  [  1/2*  ]  [  VM]/  dq^,  |  | 2  W(k,k’) 

In  the  discussion  by  Das  Sarma,  et  al  (1988),  the  overlap  integral  is  ignored. 

To  compute  equation  (2.44)  we  write  the  three  dimensional  expression  for  W(k,k’)  in 
terms  of  the  in-plane  and  perpendicular  components  q^,  we  also  recognize  that  only  the 
in-plane  phonon  vectors  take  place  in  the  selection  rule.  We  find: 

(2.45)  W(k,k’) 

=  [4x  2  ft  Fo/(mV[qy2  +  Q2  )]{Nq5  [E(k+Q)-E(k)-ftu0]  +  (Nq  + 1)5  [E(k-Q)-E(k)  +  ftu0]} 
where  the  upper  case  Q  denotes  the  in-plane  phonons: 

Ignoring  any  overlap  integral  contributions  the  integration  over  all  ^  is  direct  (assume 
contributions  from  negative  as  well  as  positive  values)  and  equal  to 

(2.46)  w(k,k’) 

=  [2*  2  ft  Fo/(nvlQ)]{Nq5  [E(k + Q)-E(k)-ftu0]  +  (Nq  + 1)5  [E(k-Q)-E(k)  +  ft«G] } 
where  w(k,k’)  retains  the  dimension  of  l/time\ 

At  this  point  we  are  interested  in  evaluating 

(2.47)  k£ft(k’-k)W(k’,k)=  krXft(k’-k)w(k\k) 

where  the  sum  on  the  right  hand  side  is  over  the  two  dimensional  quantum  states.  Thus 
we  are  interested  in  evaulating: 

(2.48)  kl£ft  (k’-k)w(k’,k)  =  [AJ{  4tt  2  )]/d 2  kft  (k’-k)w(k\k) 

As  in  the  case  of  the  three  dimensional  phonon  distribution  and  three  dimensional  ’k’ 
space,  here  we  are  going  to  assume  that  k  represents  the  polar  axis  and  that  the  phonon 
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vector  in  the  x-y  plane  submits  to  the  following  decomposition: 

(2.49)  Q  =  Q(kxcos$  +/t2sind ) 

where  are  unit  vectors.  Then,  for  phonon  absorption : 

(2.50)  S  [E(k+ Q)  -  E(k)  -  hu>Q]  =  [2mM 2  ]&  [Q2  +  2kQ  cose  ] 

or  with  ’g’  representing  the  argument  of  the  delta  function  in  equation  (2.50) 

(2.51)  6 [E(k+Q) - E(k) - ftw0]  -  [2m/ft2]S[g] 

The  zeroes  of  g(x)  are 

(2.52)  012  =  -k[cos0  ±  (cos2 5  +  huJEf1] 

where  E  =  E(k)  =  h  2  k2/2m.  In  view  of  the  constraint  that  Q  >  0,  the  only  possible  root  of 
g(x)  is 

(2.53)  Q2  =  -k[cos0  -  (cos20  +  fiuJEY'],  0£  6  £  2n 

where  we  notice  that  the  range  of  integration  is  now  over  2it,  as  is  required  by  the  two 
dimensional  nature  of  the  transport.  Thus  for  phonon  absorption  the  conservation 
condition  is: 

(2.54)  S  [E(k + Q)-E(k)-«  w  0]  =  [m/fi 2  ]S  [Q-Q  2  ]/[k(cos  2  0  +ftWo/E)X], 

0  £  e  <.  2-k 

Performing  the  integration  of  equation  (2.48)  for  absorption: 

(2.55)  kl»);*(k,*k)w(k’,k)  =  -[Kx2FoNq]0;*/2d<?[cos2tf/(cos2<?  +ftw0/Ef‘] 

We  defer  a  discussion  of  the  above  until  after  a  description  of  scattering  due  to  phonon 
emission.  The  first  thing  we  need  are  selection  rules  for  phonon  emission.  From  the 
selection  rules  for  three  dimensional  transport,  we  find: 

(2.56)  6  [E(k-Q)-E(k)  +  h  w0]  =  [2m !fi 2  ]S  [Q2  -2kQcos0  +  2mu0/ft  ] 
again  defining  a  quantity  ’g’  as  the  argument  of  the  delta  function  for  phonon  emission: 
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(2.57) 


S  [E(k-Q)-E(k)  +  ft wQ]  -  [2m /ft 2  ]S  [g] 


The  zeros  of  g 

Q3>4.=  k[cos5  ±  (cos 2B  -hij)cJEY,\ 
for  phonon  emission: 

(2.58) 

S  (E(k-Q)-E(k)  +  ft  »0]  =  [m/ft  2  }{&  [Q-Qs  ]  +  S  [Q-Q4  ]}/[k(cos2  9  -hw^Y  1 
For  E  <  ft (j)0,  there  is  no  solution  for  03  4.  For  E  >  ftu0,  and  defining 

(2.59)  =  cos’1  (ftWo/E)^ 

solutions  exist  for  in  the  first  and  fourth  quadrants,  subject  to  the  constraint 

(2.60)  cos20  >  ftujj/ E 

Thus  performing  the  integration  of  equation  (2.48)  for  phonon  emission: 

(2.61) 

kJ£ft  (k’-k)w(k’,k)  =  -[2kxF0(  Nq  + 1)][0  J*  E’£  d*  [cos2  e  /(cos2  8  -ftWo/ Ef‘] 

Combining  equations  (2.55)  and  (2.61)  the  net  scattering  rate  for  two  dimensional 
electrons  scattered  by  three  dimensional  phonons  is: 

(2.62)  r  =  [2F0/(2mE)^  ]  { Nq]  0  J*/2  dd  [cos2  5 /(cos  2  8  +  ftWo/E)X] 

(Nq  +  l)][0/tf  E'f  &  [cos2  8 /(cos2  8  -ftWo/E)^]} 

In  order  to  place  the  above  results  in  perspective  with  respect  to  the  scattering  of  three 
dimensional  electrons  by  three  dimensional  phonons  we  rewrite  equation  (238)  for  the 
rate  of  scattering  three  dimensional  electrons  by  three  dimensional  phonons: 

(2.38)  r(E)  =  F„/[2mE]X[Nq{x  +.(#Uo/E)sinh-‘(E/*Uo)K} 

+  <Nq+l){*-+  (ftUo/EJsinh-1  (E/*u0-l)8 }] 


-14- 


The  similarity  is  apparent. 


For  incorporation  of  the  above  results  into  expressions  for  evaluating  the  mobility  of  the 
carriers  equation  (2.38)  for  three  dimensional  electrons  is  inserted  into  equation  (2.34), 
whereas  equation  (2.62)  is  inserted  into  equation  (2.37).  We  have  carried  the  analysis  to 
the  point  of  integration  for  the  case  of  two  dimensional  electrons,  but  not  further.  For 
comparison  we  also  offer  the  expressions  for  the  three  dimensional  electrons.  Consider 
first  three  dimensional  electrons  and  define  xe  =  huQJ kBTe.  Then: 

(2.63)  r  (x)  =  Fo/^mkgTgx]^  [Nq{*  +  -(xo/xjsinh' 1  (x/x^  } 

+  (Nq  +  l){x'  +  (x0/x)sinh-1  (x/xq-1)*  }] 
-F0/[2mkBTex^Q(x) 

Thus  the  rate  of  momentum  loss  as  obtained  from  the  above  description  for  three 
dimensional  electrons,  when  combined  with  equation  (2.34a)  yields: 

(2.64) 


[a  njfi  kj/a  t]coli  =  -njfi  kd[4F0/3(2it  mkBTe)*  ]exp-[X  2  k^ 2  ]  0  J"  dxxQ(x)U(x)exp-x 

which  is  the  expression  used  for  computing  the  scattering  integral  for  polar  phonons  in  the 
SRA  code. 

The  comparable  expression  for  the  scattering  of  two  dimensional  electrons  by  three 
dimensional  phonons  is  obtained  with  the  scattering  rate  rewritten  as: 

(2.65)  r  =[2F0/(2mkBTex)^]{Nq]J7t/2da[cos2tf/(cos2a  -fx^x)*] 

(Nq  + 1)][  J*  e-£  da  [cos2  a /(cos 2  a  -Xg/x)*]} 

-[2F0/(2mkBTex)^]R(x) 


Then  from  equation  (2.37)  we  find: 

(2.66) 

[a  ni*  kj/a  t]colj  =  -i0nft  k(j[2F0/(2mkBTe)^]exp-[X 2  kj2  ]0  /°°dxx^R(x)X(x)exp-x 
Equation  (2.66)  is  presently  being  placed  into  a  form  suitable  for  the  MBTE  algorithm. 
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The  Density  Matrix 


The  space  charge  profiles  for  the  pseudomorphic  HEMT  were  computed  from  the 
moment  equations.  But  as  a  check  on  their  significance  as  well  as  a  means  of  addressing 
contact  problems,  as  discussed  in  later  sections,  the  liouville  equation  was  solved.  The 
Liouville  equation  is  a  quantum  mechanical  operator  equation.  In  commutator  form,  the 
equation  is: 

(2.67)  ihd  p  op/3 1  =  [H,p  0p] 

where  p0p  is  the  density  operator  and  H  is  the  system  Hamiltonian.  The  Liouville 
equation  is  expressed  with  respect  to  a  specific  representation,  the  most  popular  being  the 
mixed  representation  discussed  by  Wigner  (1932).  In  the  mixed  representation  the 
distribution  function  is  the  Wigner  function  and  quantum  transport  is  described  by  a 
differential  equation  that  is  ostensibly  similar  to  Boltzmann  equation  and  may  be  obtained 
from  the  latter  by  replacing  the  force  term  in  the  Boltzmann  equation  with  (in  one  space 
and  one  momentum  direction): 

+  ( 1/ift )( l/2x  ft  )_,  J  +  00  dp^  J  +  "  dx’f (p’pc’)[  V  (x + x’/2,t)- V  (x-x72,t)]exp[i(p-p’)x7  h  ] 

An  alternative  representation  is  that  implemented  at  SRA.  The  alternative  representation 
is  generally  referred  to  as  the  coordinate  representation;  and  the  quantity  computed  is 
called  the  density  matrix,  which  is  related  to  the  mixed  representation  (Wigner  function) 
through  (in  one  space  and  one  momentun  direction): 

(2.68)  p(xpc\t)  =  (l/(2xh)]a8/“dpf(x,p,t)exp[i(x-x’)p/ft  ] 

In  the  AFOSR  sponsored  SBIR  study  the  density  matrix  equation  of  motion  in  the 
coordinate  representation  was  solved  for  electrons. 

The  Liouville  equation  in  the  coordinate  representation  for  electrons  (ignoring  any 
spatially  dependent  effective  mass  and  dissipation)  is  the  following  differential  equation 
for  electrons  pe(xpc’,t): 

(2.69)  in  3  PfJd  t  =  -(*  2  /2m)(v 2  -v’2  )pe  +  [Ec(x,t)-Ec(x\t)]?e 

The  equation  of  motion  of  the  density  matrix  is  a  Schrodinger-like  equation.  In  its  simplest 
form  the  density  matrix  for  electrons  can  be  obtained  from  a  solution  to  Schrodingers 
equation  through  the  relation: 
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(2.70) 


p(xtx’,t)"Z/>ni1m^m(x’0^  m(x'’0 

In  equation  (2.69)  Ec(x,t)  is  the  conduction  band  energy  that  would  appear  in  the  single 
particle  Schrodinger  equation.  Ec(x,t)  is  obtained  from  the  vacuum  potential  energy 
V(x,t)  from,  the  Anderson  rule: 

(2.71)  Ec(x,t)  =  V(x,t)-X(x) 

where  X(x,t)  is  the  electron  affinity  of  the  material.  pe(x,x’,t)  is  the  density  matrix  for 
electrons  in  the  coordinate  representation. 

All  observables  such  as  density,  current,  energy,  are  obtained  from  the  diagonal 
components  of  the  density  matrix  and  suitable  derivatives  thereof  Grubin,  et  al  (1992). 
Thus  the  electron  density  is  pe(xpc): 

(2.72)  Mx)-Pe(x*x) 

The  potential  energy  of  the  vacuum  is  governed  by  Poisson’s  equation: 

(2.73)  vx[<  (x)Wx)  -  -e’  [(*,e(x)-fh(x))-(Nd  +  (x)-Pa-(x))] 

The  background  density  in  the  above  equation  is  assumed  as  a  jellium  doping  distribution 
with  N^j +  representing  the  ionized  donors  and  Pa'  representing  the  ionized  acceptors. 

Reduction  of  Liouville  Equation  for  One  Dimensional  Spatial  Transport 

In  the  calculations  described  in  this  study  we  assumed  Boltzmann  statistics,  spatial 
variations  only  along  the  x-direction,  and  free  particle  behavior  along  the  y-and 
z-directions.  The  equation  of  motion  of  the  density  matrix  was  transformed  to  center  of 
mass  ’r’  and  nonlocal  ’f  ’  coordinates: 

(2.74)  r  =  (x  +  x’)/2, 
f  =  (x— x’)/2, 

In  these  terms  the  density  matrix  is  re-expressed  as: 
p=»p(r  +  f  ,r-f ), 

and  equation  (2.69)  for  electrons  becomes: 
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(2.75) 


\h  a  p  q/8  t  =  -{h  2  /2m)d  2  pe/drd{  +  [Ec(r  +  f  ,t)-Ec(r-f  ,t)]p  e 


In  these  variables  Poissons  equation  reads: 

(2.76)  .  3/3r[<(r)aV/ar]  =  -eM(oe(r)-Oh(f))-(Nd  +  (r)-Pa-(r))) 

The  diagonal  components  of  solutions  to  equation  (2.75)  (along  the  diagonal  r=x  and 
f  =  0)  provide  the  density  for  the  electrons  with  a  similar  description  for  the  holes. 

Note,  time  independent  steady  state  conditions,  implies:  a  pjd  t  =  0 
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3.  Steady  State  Simulations  of  a  Pseudomorphic  HEMT 


Effects  of  Metalization 

The  use  of,  AlGaAs/GaAs  HEMTS  results  in  higher  electron  velocities  and  mobilities  at 
low  lattice  temperatures.  High  aluminum  mole  fractions,  however,  also  have  the 
deleterious  property  of  increased  bulk-related  trapping  centers.  To  reduce  these  traps  the 
device  designer  has  several  options;  one  is  to  reduce  the  A1  concentration,  which  has  the 
consequent  effect  of  reducing  band  offset,  and  lowering  the  2DEG  at  the  interface.  An 
alternative  was  demonstrated  by  Ketterson  et  al  (1987);  where,  for  sufficiently  thin 
InGaAs  layers,  the  compressive  elastic  strain  is  accommodated  by  tetragonal  distortion  of 
the  InGaAs  lattice  rather  than  the  generation  of  misfit  dislocations.  The  result  is  a  HEMT 
with  satisfactory  off-set  voltage  and  reduced  bulk-related  trap  centers.  HEMTs 
constructed  from  such  a  quantum  well  are  generally  called  pseudomorphic  HEMTs. 

The  fabrication  of  a  HEMT  generally  involves  the  placement  of  contacts  on  a  narrow 
band-gap  material  such  as  gallium  arsenide  (or  indium  gallium  arsenide),  which  in  turn  is 
deposited  on  the  wide  band  gap  material.  Simulating  such  a  structure  is  always  difficult  as 
it  requires  dealing  with  transport  through  a  barrier,  necessitating  the  implementation  of 
quantum  transport  algorithms.  With  regard  to  the  issue  of  transport  through  the  barrier  it 
is  recognized  that  there  is  a  body  of  opinion  suggesting  that  metalization  destroys  the 
integrity  of  the  barrier.  Rather  than  ignore  the  barrier  issue  we  asked.  What  if?,  the  barrier 
is  present. 

This  issue  was  examined  generically,  through  implementation  of  quantum  algorithms 
developed  in  part  with  support  from  AFOSR.  If  we  assume,  for  example  an  FET  structure 
whose  topmost  portion  can  be  represented  by  figure  3.1,  then  it  is  possible  to  implement 
the  SRA  quantum  transport  algorithm  and  ask  what  if?  questions  concerning  the  transport 
of  carriers  to  the  quantum  well  for  the  case  when  the  integrity  of  the  heterostructure  layers 


Narrow  Band  Gap  Material 


Wide  Band  Gap  Material 


Narrow  Band  Gap  Quantum  Well 


Fig  3.1a  Generic  top  layers  of  a  HEMT  structure. 
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is  not  compromised  by  the  contact  metalization  process. 

The  SRA  quantum  transport  simulations  involve  implementing  an  algorithm  for  solving 
the  Liouville  equation  in  the  coordination  representation.  Figure  3.1b  displays  the 
self-consistent  potential  and  charge  distributions  (solutions  to  Poisson’s  equation  are 
included),  under  equilibrium  conditions,  for  a  uniformly  doped,  350A  wide  band  (barrier 
height  of  300  mev),  region  surrounded  by  two  narrow  band  gap  regions.  Such  regions  may 
be  expected  to  mimic  under  equilibrium  the  source  and  drain  regions  of  a  pseudomorphic 
HEMT  with  uniform  doping  through  the  quantum  well. 

Several  features  are  highlighted:  (i)  a  two  dimensional  electron  gas  forms  on  both  sides  of 
the  wide  band  gap  region;  (ii)  the  barrier  as  seen  by  carriers  entering  the  structure  is 
reduced  by  nearly  10  mev,  and  is  replaced  by  two  small  tunneling  regions  at  the  ends  of  the 
barrier.  The  result  clearly  suggests  that  a  significant  component  of  current  from  the 
source  contact  to  the  quantum  well  can  be  accomplished  by  carriers  tunneling  through  the 
thin  barriers,  as  well  as  by  thermionic  emission  over  the  barrier.  Note,  the  considerably 
reduced  free  carrier  concentration  in  the  barrier  region.  On  the  basis  of  this  calculation 
there  is  likely  to  be  an  insignificant  amount  of  parasitic  source  to  drain  current  flowing  in 
the  wide  bandgap  material  parellel  to  the  2DEG. 


-GOO  -300  0  300  600 

Distance  (Angstroms)  Along  The  Diagonal 


-600  -  300  0  300  600 

Distance  (Angstroms)  Along  The  Diagonal 


Figure  3.1b.  Self  consistent  distribution  of  potential  energy  and  charge  for  a  uniformly  doped 
narrow  band  gap-wide  band  gap-  narrow  band  gap  structure. 


The  situation  when  the  doping  in  the  structure  ends  within  the  wide  gap  material  is 
qualitatively  similar  on  the  contact  side  of  the  structure.  However,  at  the  interface 
between  the  wide  bandgap  region  and  the  quantum  well  the  density  of  the  2DEG  is 
reduced,  as  shown  in  figure  3.1c 
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Figure  3.1c.  As  in  figure  3.1b,  but  where  the  doping  ends  at  the  quantum  well. 


The  Phase  II  calculations  that  were  performed  for  the  PHEMT  were  based  upon 
semiclassical  physics,  and  tunneling  contributions  were  not  included.  A  different  study  on 
HEMTs  which  included  semiclassical  physics  permitted  thermionic  emission  over  the 
barrier  and  indicated  that  the  barrier  provided  an  unrealistically  high  source  resistance. 
This  necessitated  performing  calculations  in  which  boundary  conditions  were  imposed 
directly  on  the  wide  band  gap  material,  a  feature  common  to  most  HEMT  studies. 
Because  of  the  availability  of  the  quantum  transport  algorithm  at  SRA,  the  potential  and 
mobile  charge  density  distribution  were  computed  for  a  situation  in  which  boundary 
conditions  were  imposed  directly  on  the  wide  band  gap  material,  as  was  done  in  the  full 
PHEMT  simulations.  These  calculations  provided  some  measure  of  the  nature  of  the 
semiclassical  approximations  we  made. 


Figure  3.2a  displays  a  sequence  of  calculations  for  a  1200  A  structure.  The  left  600  A  is 
comprised  of  a  wide  band  gap  material  (300  mev  barrier);  the  right  side  is  a  narrow 
bandgap  material.  The  sequence  of  calculations  is  for  a  doping  distribution  of 
101  */101  a,  101  V101 7,  101  */101 6 ,  101  VIO1  8,  respectively.  There  are  several  points  to 
note:  (i)  approximately  150  A  of  the  wide  band  gap  material  is  depleted  of  carriers;  (ii) 
the  2DEG  forms  within  the  first  100  A  of  the  narrow  bandgap  material;  (iii)  the  density  of 
the  2DEG  is  insensitive  to  quantum  well  doping  levels  below  101  6 /cm3.  The  potential 
distribution  is  also  shown,  and  incorporates  the  conduction  band  offset.  Note  that 
approximately  2/3  of  the  offset  voltage  falls  across  the  wideband  gap  material.  How  do 
these  calculations  compare  to  the  semiclassical  calculations? 


Semiclassical  calculations  performed  under  the  Phase  I  study,  in  which  metallization  is 
assumed  computationally  to  occur  at  the  wide  band  gap  material  show  the  distribution 
displayed  in  figure  3.2b.  In  comparing  figures  3.2a  and  3.2b,  the  first  point  to  note  is  that 
the  density  is  higher  in  the  classical  calculation  (in  the  quantum  mechanical  calculation 
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Figure  3.2a  Quantum  mechanical  calculations  of  the  distribution  of  charge  and  potential  for 
a  wide  band  gap!  narrow  band  gap  structure  with  varying  doping  concentrations. 


x(am) 

Figure  3.2b.  Distribution  in  charge  under  the  source  region  for  a  400  A  wide  band  gap  region 
adjacent  to  a  narrow  band  gap  region.  From  Phase  /  Final  Report,  figure  21  a. 


variations  in  the  effective  mass  were  not  taken  into  account).  The  second  point  to  note  is 
that  the  density  drops  precipitously  at  the  end  of  the  quantum  well,  a  feature  not 
accounted  for  in  the  quantum  calculations  displayed  here,  although  a  rapid  decrease  in 
charge  at  the  end  of  quantum  well  has  been  observed  by  the  investigators  in  other 
calculations.  The  broad  conclusion  that  can  be  drawn  from  a  comparison  of  figures  3.2a 
and  3.2b  is  that  the  parasitic  current  flow  parallel  to  the  2DEG  but  in  the  wide  band  gap 
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material  is  likely  to  be  overestimated  in  the  semiclassical  calculations. 

The  Structure  Simulated 

The  structure  that  was  simulated  under  the  Phase  II  study  is  shown  in  figure  3.3.  Note  that 
boundary  conditions  were  imposed  on  the  wide  band  gap  material.  The  simulations  were 
carried  out  for  a  range  of  bias  conditions  using  the  moments  of  the  Boltzmann  transport 
equation  algorithm.  The  device  consists  of  0.5  pm  source  and  drain  contacts  with  a  source 
drain  spacing  of  1.0  pm.  An  0.25  /im  gate  is  placed  symmetrically  between  the  source  and 
drain.  The  device  structure  consists  of  a  350A  layer  of  AlGaAs  doped  to  1  x  lO^/cm^, 
followed  by  a  30A  undoped  layer  of  AlGaAs.  Below  this  is  a  150A  deep  channel  of 
undoped  InGaAs,  which  sits  on  a  GaAs  substrate. 

The  Current  Voltage  Characteristics 

The  Phase  I  current  voltage  characteristics  for  this  structure  with  a  device  width  of  fifty 
microns  obtained  through  solution  to  drift  and  diffusion  equations  are  shown  in  figure  3.4, 
while  those  obtained  using  the  moments  of  the  Boltzmann  transport  equation  are  shown  in 
figure  3.5.  For  the  MBTE  calculation,  gate  bias  levels  of  0.4,  0.0,  and  -1.0  volts  were 
considered  over  a  drain  bias  level  between  0.0  and  3.0  volts.  The  current  voltage 
characteristics  of  the  device,  as  predicted  by  the  *  "ITE  simulation  code  show  current 
levels  that  are  of  the  order  of  three  times  la.ger  than  that  obtained  with  the  dirft  and 
diffusion  results,  a  feature  that  is  consistent  with  effects  of  velocity  overshoot.  There  is, 
however,  structure  in  these  curves  that  is  not  present  in  the  drift  and  diffusion  calculations. 

Internal  Dynamics 

The  current  voltage  characteristics  represent  a  coarse  summary  of  the  internal  dynamics 
of  the  PHEMT.  It  is  the  internal  dynamics  that  is  discussed  next  through  line  and  contour 
plots  of  the  internal  distribution  of  carriers,  temperature  and  potential.  Transient  effects 
associated  with  switching  between  bias  points  is  also  discussed  in  terms  of  these  plots. 

Interpretation  of  these  results  is  aided  by  examining  Figure  3.6,  which  shows  a  portion  of 
the  computational  mesh  used  in  the  simulation.  The  location  of  the  material  and  doping 
interfaces  are  indicated  in  the  figure.  The  domain  shown  includes  only  approximately  the 
first  425A  of  the  substrate.  As  will  become  evident  in  the  subsequent  figures,  the 
remainder  of  the  substrate  plays  no  significant  role  in  the  device  dynamics.  The  points 
labeled  1  through  3  represent  selected  points  where  the  transient  history  of  a  variety  of 
quantities  will  be  examined. 
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Figure  3.3.  Schematic  of  the  Pseudomorphic  HEMT  (Device  Width  =  50  microns). 
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12 


Figure  3.4.  Current  voltage  characteristics  for  a  50  micron  wide  PHEMT  as 
obtained  from  the  drift  and  diffusion  equations.  From  Phase  I  Final  Report 


Figure  3.5  Current  voltage  characteristics  for  a  50  micron  wide  PHEMT  as  obtained  from 
the  moments  of  the  Boltzmann  transport  equation. 
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Figure  3.6  Computational  grid  for  the  PHEMT. 


Equilibrium  Calculation 

Figure  3.7  shows  the  potential  distribution  (a)  and  electron  density  (b)  when  V^s  and  Vgs 
are  zero.  Under  such  conditions  the  carriers  are  in  static  equilibrium  with  the  lattice  at 
300°  K.  The  contours  show  the  anticipated  symmetric  distribution  of  these  quantities. 
The  high  field  region  under  the  gate  (evidenced  by  the  closely-spaced  potential  contours) 
and  the  depletion  of  carriers  from  the  heavily-doped  AIGaAs  under  the  gate  are  clearly 
evident.  Also  evident  is  the  accumulation  (darkened  regions)  of  carriers,  a 
two-dimensional  gas,  at  the  interface  between  the  AIGaAs  and  InGaAs  materials  under 
the  source  and  drain  regions  and  their  confinement  to  the  InGaAs  channel.  The  density  in 
this  region  reaches  levels  in  excess  of  3  x  10^/cm^.  The  depletion  of  electrons  from  the 
channel  under  the  gate  is  also  evident. 

Figure  3.8  displays  line  plots  of  the  potential  and  electron  density.  The  plots  are  along  a 
line  normal  to  the  center  of  the  gate.  The  potential  plot  represents  only  the  self  consistent 
potential  as  obtained  from  Poisson’s  equation.  (The  self  consistent  potential  plus  the 
heterostruture  potential  are  incorporated  into  the  calculations  discussed  in  figures  3.2.  ) 
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Here,  Nt  represents  the  lower  valley  carriers  and  +N2  the  total  carrier  density,  the 
difference  being  upper  L-valley  electrons.  For  this  calculation  most  of  the  carriers  are  in 
the  r -valley.  This  figure  clearly  shows  the  depletion  of  carriers  directly  under  the  gate  and 
in  the  channel.  However,  the  carrier  density  in  the  channel,  due  to  accumulation  at  the 
material  interface,  remains  at  a  value  of  nearly  1.5  x  101 7 /cm3  and  then  decreases  rapidly 
in  the  substrate.  Recall  that  under  the  source  and  drain  the  density  in  the  channel  exceeds 
3X101 8 /cm3. 


b)  Electron  Density. 


Figure  3. 7  Computed  Contours  for  V^s  =  0.0  and  VqS  =  0.0 
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Figure  3.8  Computed  Distributions  Normal  to  Center  of  Gate  for  Vjs =0.0  and  V„s  = 0.0 


Nonequilibrium  Calculations 


Steady  State,  Vjs  =  1.0  v,  VgS  =  0.0  v 

The  next  set  of  figures.  Figures  3.9  through  3.12,  show  the  results  for  =  1.0,  VgS=0.0. 
In  Figure  3.9a  the  potential  distribution  shows  the  effect  of  the  applied  bias.  The 
potential  contours  under  the  gate  are  distorted  with  a  larger  potential  drop  under  the 
drain  side  of  the  gate  contact,  than  on  the  source  side  of  the  gate  contact.  There  is  also  the 
formation  of  a  longitudinal  field  from  the  source  to  the  drain  contact,  as  represented  by 
vertical  contours.  The  high  field  at  the  drain  end  of  the  gate  extends  across  the  channel  and 
into  the  substrate,  where  it  has  the  effect  of  distorting  the  charge  distribution,  without  any 
significant  influence  on  the  transport  of  carriers  to  the  contacts. 

The  density  contours,  Figure  3.9b,  when  compared  to  Figure  3.7,  show  an  extension  of  the 
depletion  region  outward  into  the  channel  and  a  reduction  of  the  carrier  density  in  the 
channel  under  the  gate.  This  occurs  because  electrons  are  drawn  from  this  region  by  the 
applied  bias.  Thus,  even  though  a  current  flows,  the  carrier  density  is  lower  in  the  channel 
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than  when  Vjs  =  0.0. 


The  current  path  from  source  to  drain  is  shown  by  the  current  streamlines  in  Fig.  3.9c. 
Here  we  observe  that  the  current  path  flows  directly  from  the  source  to  the  2  DEG.  The 
majority  of  the  current  then  flows  within  the  2  DEG  to  the  drain  side  of  the  device,  where 
it  is  collected  at  the  drain.  The  streamlines  which  intersect  the  bottom  of  the  figure 
indicate  current  paths  deeper  in  the  substrate.  These  represent  a  small  fraction  of  the 
current  flow  from  the  source  to  the  drain  contact. 

Figures  3-10  through  3-12  show  profile  plots  along  lines  normal  to  the  center  of  the 
source,  gate  and  drain,  respectively  for  the  case  V^s  =  1.0  and  Vgs  =  0.0.  For  the  source 
and  drain  region,  the  self  consistent  potential  variation  shows  the  expected  shape 
(potential  decreases  toward  the  quantum  well)  resulting  from  diffusion  of  carriers  into  the 
undoped  region  and  the  accumulation  of  carriers  at  the  AlGaAs/InGaAs  interface  under 
the  source  and  drain.  The  shape  of  the  potential  under  the  gate  corresponds  to  depletion 
under  the  gate. 

The  density  profiles  under  the  source  and  drain.  Figures  3-10b  and  3-12b,  show  the 
diffusion  of  carriers  into  the  undoped  AlGaAs  region  from  the  heavily-doped  AlGaAs 
layer  and  the  sharp  peak  in  concentration  in  the  2  DEG.  Several  points  are  noteworthy: 
(1)  The  distribution  of  charge  under  the  source,  which  may  be  expected  to  represent 
equilibrium  conditions,  differs  from  the  result  obtained  quantum  mechanically  (as 
discussed  earlier).  (2)  Under  the  influence  of  the  drain  bias  the  peak  density  in  the  2 
DEG  under  the  drain  approaches  4X101  8 /cm3 .  This  is  a  result  of  two  factors.  First,  due 
to  the  drain  bias  electrons  are  pulled  towards  the  drain  end  of  the  channel.  Second, 
electrons  easily  enter  the  2  DEG  at  the  source,  but  the  interface  between  the  AlGaAs  and 
InGaAs  acts  as  a  barrier  at  the  drain.  Electrons  pile  up  at  this  barrier  in  an  attempt  to 
reach  the  drain.  Under  the  gate,  in  Figure  3-1  lb,  the  electron  density  in  the  channel  is 
much  lower,  at  a  level  near  5xl016/cm3,  than  the  equilibrium  value  of  near 
1.5xl0x  7 /cm3. 

The  central  valley  electron  temperature  shown  in  these  three  figures  indicates  that 
electrons  enter  and  leave  the  device  nearly  at  equilibrium  with  the  lattice.  The  lattice 
temperature  is  at  the  ambient  at  the  source  and  drain  contact.  A  small  increase  in 
temperature  to  1.4  times  the  lattice  temperature,  or  420°  K  occurs  under  the  gate,  due  to 
the  high  source-drain  field  in  this  region.  This  increase  is  fairly  uniform  across  the  entire 
device,  consistent  with  the  potential  variation  shown  in  Figure  3.9a.  Within  the  context  of 
the  electron  temperature  model  used  in  the  moment  formulation,  increased  election 
temperature  implies  electron  transfer  from  the  r  valley  to  the  L  valley.  The  temperature 
rise  for  this  bias  level  was  not  great  enough  to  cause  significant  electron  transfer. 


-29- 


The  central  valley  velocity  component  parallel  to  the  heterostructure  interface  is  shown  in 
Figure  3-1  Id.  Again  the  distribution  is  fairly  uniform  across  the  device,  reaching  a  value 
of  approximately  6xl07  cm/sec.  Since  the  only  place  any  significant  numbers  of  carriers 
are  present  is  at  the  2DEG,  only  a  small  current  density  emerges. 
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b)  Electron  Density. 
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Figure  3.10  Computed  Distributions  Normal  to  Center  of  Source  for  V^s  =  1.0  and  VgS 
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Steady  State,  Vjs  =  2.0  v,  VgS  =  0.0  v 

When  the  drain  bias  is  increased  to  V(js  =  2.0  volts  while  holding  the  gate  bias  at  0.0,  the 
distribution  of  potential,  density  and  current  streamlines  is  as  shown  in  Figure  3.13.  An 
increase  in  the  distortion  of  the  potential  distribution  under  the  gate,  relative  to  the 
Vds  =  0.0  and  1.0  volt  cases,  is  observed,  as  is  the  higher  field  region  over  a  wider  area  on 
the  drain  side  of  the  gate.  Due  to  this  higher  field,  the  velocities  in  the  channel  are  greater 
and  the  carrier  density  continues  to  decrease  under  the  gate,  as  shown  in  Figure  3.13b. 
The  injection  of  carriers  into  the  substrate  is  also  increased,  as  is  indicated  by  the  density 
contours,  Figure  3.13b,  and  the  current  streamlines,  Figure  3.13c.  Note  that  a  more 
substantial  fraction  of  the  current  now  flows  through  the  substrate. 

The  distributions  under  the  source,  gate  and  drain  are  shown  for  this  bias  level  in  Figures 
3.14-3.16.  The  potential  variation.  Figures  3.14a  -  3.16a,  show  that  under  the  source  only  a 
small  field  is  present  in  the  region  of  the  AlGaAs/InGaAs  interface,  as  was  the  case  for 
Vds  =  10.  No  significant  field  appears  in  the  substrate  in  the  direction  normal  to  the 
source.  However,  under  the  gate,  Figure  3.15a  shows  that  a  field  forcing  electrons  towards 
the  substrate  exists,  and  this  field  is  responsible  for  the  current  paths  extending  further 
into  the  substrate.  Under  the  drain,  the  field  in  the  substrate  is  reversed,  drawing 
electrons  back  towards  the  channel  and  to  the  drain. 

The  density  profiles  shown  in  Figures  3.14b  -  3.16b  are  similar  to  the  V(jg  =  1.0  volt  case. 
However,  we  note  that  the  peak  in  carrier  density  under  the  source  is  reduced,  while  the 
peak  under  the  drain  increases.  This  is,  of  course,  consistent  with  the  higher  applied  drain 
bias.  Under  the  gate,  the  reduced  peak  density  is  evident,  as  is  the  longer  tail  of  the 
profile,  showing  the  injection  of  carriers  into  the  substrate.  Note  that  no  significant 
electron  transfer  has  yet  occurred.  The  central  valley  carrier  temperature  under  the  gate 
has  risen  to  1.8  times  the  lattice  temperature,  or  540®  K.  This  is  not  high  enough  to  result 
in  significant  transfer.  Finally,  we  see  that  the  central  valley  velocity  has  risen  to  about 
8xl0^cm/sec  under  the  gate. 

For  the  VgS  =  0.0  volt  cases  just  examined  it  is  apparent  that  for  of  two  volts  and 
below  no  significant  heating  of  the  carriers  occurs  and  the  velocity  in  the  channel  under 
the  gate  increases  in  such  a  manner  that  when  combined  with  the  carrier  density  profile, 
the  current  shows  a  relatively  linear  characteristic  vs  V^g.  As  V^g  increases  further,  to  3.0 
volts,  some  transfer  occurs  and  the  current  begins  to  saturate  according  to  the  simulations. 
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Figure  3.14  Computed  Distributions  Normal  to  Center  of  Source  for  V^s  =  2.0  and  Vgs 
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Steady  State,  V^s  =  2.0  v,  VgS  =  -1.0  v 

When  the  gate  bias  was  reduced  to  -1.0  volts,  the  drain  current  assumes  negligible  values 
for  V^s  below  2.0  volts.  Contour  plots  of  the  potential  and  electron  density  for  this  case 
are  shown  in  Fig.  3.17.  As  shown  in  Figure  3.17a  the  potential  and  carrier  distributions 
under  the  gate  show  an  extended  depletion  region  penetrating  deep  into  the  substrate. 
The  channel  is  void  of  carriers  under  the  gate  and  the  device  is  effectively  pinched  off. 
This  result  is  re-enforced  in  Figure  3.18,  where  profile  plots  under  the  gate  are  shown. 
The  lack  of  any  free  electrons,  even  2,000A  into  the  substrate,  clearly  shows  there  will  be 
negligible  current.  Due  to  the  applied  drain  potential,  the  velocity  of  what  few  electrons 
do  flow  is  at  a  level  of  2xl0^cm/sec,  but  the  net  flux  of  electrons  is  approximately  zero. 

StecAy  State,  Vjs  =  1.0  v,  VgS  =  0.4  v 

When  the  gate  bias  is  raised  to  0.4  volts  and  the  drain  bias  set  at  1.0  volts,  the  contours  of 
potential,  density  and  current  streamlines  appear  as  shown  in  Figure  3.19.  These  results 
may  be  compared  to  Figure  3.9,  for  Vgs  =  0.0  volts.  The  potential  contours  are  more 
widely  separated  and  the  source  to  drain  field  increases  more  gradually  than  for  the 
Vgs  =  0.0  v  case.  It  is  also  apparent  that  there  is  a  greater  potential  drop  between  the 
InAlGa  side  of  the  channel  and  the  drain  contact.  Since  a  higher  current  flows,  more 
electrons  must  cross  this  barrier  to  reach  the  drain  and  more  energy  is  required.  The 
density  contours  clearly  show  the  reduced  depletion  under  the  gate  contact  region.  It  is 
also  evident  from  the  extent  of  the  heavily-shaded  areas  in  this  figure  that  the  electron 
concentration  in  the  2  DEG  is  much  greater  at  the  drain  end  of  the  device.  This  is  similar 
to  the  Vgs -0.0  volt  case,  but  more  extreme.  The  current  streamlines  while  qualitatively 
similar  to  the  VgS  =  0  case,  showing  that  the  current  flows  primarily  in  the  2  DEG  region, 
also  show  a  significant  amount  of  current  crowding  on  the  gate  side  of  the  drain  contact, 
leading  to  the  possibility  of  premature  breakdown  in  this  region. 

Of  the  profile  plots  normal  to  the  gate  shown  in  Figure  3.20,  it  is  of  particular  interest  to 
examine  the  electron  density.  With  0.4  volts  forward  bias  on  the  gate,  the  density  in  the  2 
DEG  approaches  2.5x10 1 7 /cm 3  and  it  is  apparent  that  the  electrons  in  the  heavily-doped 
AlGaAs  layer  are  not  fully  depleted,  as  in  the  cases  when  VgS  =  0.0  v  or  below.  Also  note 
that  the  electron  temperature  under  the  gate  is  much  higher  than  at  vgs  =  0.0v.  Higher 
drain  voltage  results  in  further  increases  in  electron  temperature,  and  consequent  electron 
transfer,  and  appears  to  be  the  largest  contributor  to  the  slope  change  in  the  current 
voltage  characteristics,  as  displayed  in  figure  3.5. 
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Figure  3.18  Computed  Distributions  Normal  to  Center  of  Gate  for  Vjs  =  1.0  and  VgS  -  - 1.0 
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d)  Gamma  Valley  Horizontal  Velocity  Magnitude. 


Figure  3.20  Computed  Distributions  Normal  to  Center  of  Gate  for  V^s  - 1.0  and  VgS  ~  0. 


Steady  State,  V^s=3.0  v,  VgS  =  0.4  v 

When  the  drain  voltage  is  raised  to  3.0  volts,  with  Vgs  =  0.4,  the  results  are  as  shown  in 
Figure  3.21.  The  effect  of  the  increased  drain  bias  dominates  the  potential  distribution  in 
the  channel,  creating  a  broad  region  of  high  field  in  the  source  drain  direction.  The 
density  distribution  shows  a  reduction  in  the  density  levels  at  the  source  end  of  the  2  DEG 
and  under  the  gate.  Additionally,  at  the  drain  side  of  the  device  the  electrons  accumulate 
at  the  extreme  right  side  of  the  device,  under  the  drain.  The  current  streamlines  show  a 
deeper  penetration  into  the  substrate,  which  is  characteristic  of  increasing  the  drain  bias 
level,  and  further  current  crowding. 

Distribution  normal  to  the  gate,  shown  in  Figure  3.22,  indicates  substantial  electron 
transfer  and  a  central  valley  electron  temperature  near  1260°  K.  Further  there  is  a  large 
spread  of  the  charge  into  the  quantum  well.  While  there  is  a  reduction  in  the  net  charge 
under  the  gate  the  central  valley  velocity  reaches  8xl07  cm/sec,  and  is  in  major  part 
responsible  for  the  increased  current  flow  at  this  bias  level.  The  high  degree  of  electron 
transfer,  due  to  the  high  current  flow  and  increased  drain  bias,  gives  rise  to  the  current 
saturation  observed  in  the  current  voltage  characteristics  shown  in  Figure  3.5.  In  addition, 
it  is  pointed  out  that  when  electron  transfer  occurs  there  is  a  substantial  increase  in  the 
effective  mass.  Quantum  transport  calculations  indicate  that  this  increase  in  effective 
mass  reduces  the  population  of  the  quantized  two  dimensional  electron  gas.  While  this 
quantization  is  not  incorporated  in  this  study;  a  reduction  in  the  current  levels  should 
accompany  this  change. 
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c)  Current  Streamlines. 


Figure  3.21  Computed  Contours  for  V^s  -3.0  and  VgS  =  0.4 


c)  Gamma  Valley  Temperature. 


d)  Gamma  Valley  Horizontal  Velocity  Magnitude. 


Figure  3.22  Computed  Distributions  Normal  to  Center  of  Gate  for  Vjs =3.0  and  VgS  =  0.4 


4.  Transient  Simulations  of  the  Pseudomorphic  HEMT 

Introduction 

Perhaps  one  of  the  most  significant  features  of  the  ability  to  perform  transient  calculations 
with  either  the  moments  of  the  Boltzmann  transport  equation  or  Monte  Carlo  methods  is 
that  the  effects  of  velocity  overshoot  on  the  transient  behavior  of  the  devices  can  be 
determined.  It  is  recalled  that  in  the  early  seventies  the  initial  work  of  Ruch  (1972) 
suggested  the  possibility  of  very  short  time  scale  devices.  Subsequent  work  by  others,  (see 
e.g.,  Grubin,  et  al,  1982)  indicated  that  high  speeds  may  not  necessarily  be  attainable  since 
the  mechanism  for  relaxation  involved  several  time  scales,  (i)  the  momentum  relaxation 
time  and  (ii)  the  energy  relaxation  time.  There  is  also  an  intervalley  intervally  scattering 
rate  to  contend  with.  Finally  there  is  the  fact  that  it  takes  a  finite  amount  of  time  to  charge 
a  capacitor  and  that  the  displacement  currents  must  be  accounted  form.  In  an  earlier 
study,  Grubin  and  Kreskovsky  (1985)  demonstrated  that  on  the  time  scales  of  interest  the 
transients  associated  with  the  displacement  currents  are  of  the  same  order  as  the  transient 
hot  electron  time  scales.  Thus  the  simplest  question  to  ask  is  :  How  will  the  PHEMT 
respond  to  sudden  and/or  controlled  changes  in  the  time  dependent  applied  bias?  This 
question  is  addressed  below. 

The  first  set  of  calculations  discussed  involves  the  response  of  the  PHEMT  structure  to 
sudden  changes  in  bias.  This  is  followed  by  a  set  of  calculations  involving  a  controlled  rate 
of  change  of  bias. 

VgS  =  +  0.4  v;  Drain  Bias  is  Switched  from  1.0  vtoZOv  in  Zero  Time 

Several  transient  simulations  were  performed  to  assess  the  effects  of  switching  of  the  drain 
bias.  The  first  calculation  involves  a  switch  in  the  drain  bias  from  1.0  volt  to  ZO  volts,  while 
the  gate  bias  was  maintained  at  0.4  volts.  The  response  of  the  contact  currents  is  shown  in 
Figure  4.1. 

Figure  4.1a  shows  the  response  of  the  total  current  (conduction  plus  displacement)  at  the 
source,  gate  and  drain.  Note  that  the  peak  in  the  drain  current  is  clipped  in  the  plot.  The 
convention  used  here  is  that  positive  gate  and  drain  current  represents  electrons  flow  out 
of  the  structure.  Negative  current  at  the  source  represents  electrons  flowing  into  the 
device.  At  1.1  psec  the  gate  current  is  approximately  zero  while  the  source  and  drain 
currents  are  approximately  equal  in  magnitude  but  opposite  in  sign. 

While  the  response  was  not  carried  out  to  a  completely  steady  result,  the  current  transients 
indicate  several  features.  The  overall  transient  has  a  time  scale  on  the  order  of  several 
picoseconds  and  appears  to  be  governed  by  the  gate  capacitance.  From  the  gate  response,  it 
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is  apparent  that  the  gate  current  displays  an  initial  transient  with  a  characteristic 
oscillation  at  a  frequency  of  the  order  of  two  THz.  The  pattern  of  current  flow  indicates, 
for  this  transient,  that  there  is  strong  current  and  charge  coupling  through  the  interaction 
of  the  gate  and  drain  regions.  However,  the  drain,  which  closely  follows  the  gate  response, 
also  sustains  a  higher  frequency  oscillation  of  the  order  of  ten  THz;  this  higher  frequency 
response  is  a  consequence  of  its  interaction  with  the  source  region. 

Examination  of  the  displacement  currents  in  Figure  20b  indicates  the  origin  of  these  high 
frequency  components  of  the  currents  are  the  rapidly  changing  levels  of  charge  in  the  2  DEG 
under  the  source  and  drain.  As  was  shown  in  previous  figures,  the  electron  distribution 
exhibits  large  peaks  at  the  material  interfaces  under  the  source  and  drain.  The  magnitude 
of  these  peaks  was  shown  to  be  highly  dependent  on  the  drain  bias  level.  Thus,  a  step 
change  in  drain  bias  results  in  a  readjustment  of  the  charge  distribution  in  this  region, 
giving  rise  to  these  capacitance  effects. 

Further  evidence  of  this  effect  is  presented  in  Figure  4.2,  where  the  time  variation  of 
various  parameters  at  a  given  point  within  the  device  are  shown.  In  Figure  4.2  the  plots 
are  labeled  "point  no.  1",  etc.  This  notation  refers  to  the  location  within  the  device  where 
the  data  was  recorded.  These  locations  are  indicated  in  Figure  3.6.  Point  1  is  within  the 
InGaAs  channel  under  the  source,  point  2  within  the  channel  near  the  drain  end  of  the  gate, 
and  point  3  within  the  channel  under  the  drain.  The  results  at  points  1  and  3,  under  the 
source  and  drain,  respectively,  show  that  there  is  very  little  transient  effect  on  the  central 
valley  temperature,  Tl,  or  velocity  magnitude,  Ql.  However,  the  central  valley  density, 
Nl,  shows  a  significant  transient  effect  with  variations  on  the  time  scale  associated  with 
the  displacement  currents  at  the  source  and  drain.  These  transients  settle  out  in  about  a 
picosecond.  Under  the  gate  a  more  significant  effect  is  observed.  The  velocity  magnitude, 
Ql,  increases  significantly  and  the  density  drops.  The  potential  at  this  point  also  rises 
approximately  0.6  volts.  Recall  that  point  2  is  located  in  the  high  field  region,  on  the  drain 
side  of  the  gate.  While  these  three  quantities  vary  dramatically,  it  is  observed  that  the 
transients  are  smooth  without  the  high  frequency  components  observed  at  the  source  and 
drain.  This  is  because  we  are  dealing  with  transients  associated  with  length  scales  along 
the  channel,  i.e.,  the  source-drain  spacing.  Under  the  source  and  drain,  the  relevant 
length  scale  is  that  associated  with  the  accumulation  layer  at  the  material  interface,  i.e., 
the  2  DEG  channel  depth.  We  note  also  that  at  point  2  the  transients  also  decay  in 
approximately  a  picosecond,  with  the  exception  of  the  central  valley  temperature.  As  we 
shall  see,  the  time  scale  for  the  temperature  to  adjust  to  the  bias  change  is  somewhat 
longer,  on  the  order  of  3  picoseconds.  It  is  this  3  picosecond  time  scale  that  determines  the 
time  to  relaxation  for  the  PHEMT. 
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Figure  4.1  T ransient  Response  for  a  Change  in  Bias,  V^s  =  1.0  to  2.0  and  V. 
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Figure  4.2  Transient  Response  for  a  Change  in  Bias,  Vfo  =  1.0  to  2.0  and  V„s  ~  0.4 


VgS  =  +  0.4  v;  Drain  Bias  is  Switched  from  1.5  v  to  3.0  v  in  Zero  Time 

The  results  of  a  similar  transient  simulation,  but  with  a  greater  change  in  the  drain  voltage, 
from  1.5  to  3.0  volts,  are  shown  in  Figures  4.3  and  4.4.  These  figures  exhibit  the  same 
characteristics  and  time  scales  discussed  previously.  However,  this  simulation  was 
extended  in  time  to  3.0  picoseconds.  In  Figure  4.4  the  temperature  variation  at  point  2 
clearly  reaches  a  new  steady  value  after  3.0  picoseconds.  The  results  shown  in  Figures 
4. 1-4.4  show  that  switching  time  of  the  device  is  controlled  by  the  thermal  energy 
relaxation  time.  This  is  not  unexpected,  since  these  relaxation  times  are  the  longest. 

What  is  the  significance  of  the  above  result?  It  is  recalled  that  during  the  Phase  I  study  it 
was  argued  that  the  speed  of  a  device  is  the  time  it  takes  the  device  response  to  fall  within 
5%-10%  of  its  steady  state  response  corresponding  to  those  bias  conditions.  The  results  of 
the  Phase  II  study  indicate  that  for  this  specific  structure,  the  time  determining  the 
response  is  the  longest  response  time,  which  is  the  energy  relaxation  time.  Since 
momentum  relaxation  occurs  on  a  much  shorter  time  scale  than  the  relaxation  of  the 
electron  temperature,  energy  relaxation  in  this  study  is  dominated  by  temperature 
relaxation.  In  the  above  calculations,  energy  relaxation  is  represented  as  being  somewhat 
greater  than  3  ps.  As  a  consequence  of  this  if  a  signal  is  placed  on  one  of  the  contacts  with 
a  period  that  is  of  the  order  of  the  energy  relaxation  time  or  shorter,  that  the  transport 
quantities  will  initially  deviate  from  their  steady  state  values,  and  ultimately  will  adjust  to 
periodic  values  that  will  differ  signficantly  from  their  steady  state  values.  Figure  4.5 
displays  the  results  of  such  a  calculation. 

Vds  =  3.0  v,  VgS(t)  =  VgS( t  +  5. Ops);  VgS  =  0.4 ±  0.4[l-exp-t/(lps)],  0  <t  <  2.5ps. 

Figure  4.5  displays  the  response  of  the  PHEMT  to  a  exponentially  rising  and  decaying  gate 
voltage.  Figure  4.5a  is  a  sketch  of  the  applied  gate  bias  which  is  periodic  with  a  period  of 
the  order  of  the  energy  relaxation  time.  The  total  current  through  the  three  contacts  is 
displayed  in  figure  4.5b.  Note  that  the  source  current  appears  to  saturate  after  1.0  ps, 
while  the  gate  voltage  continues  to  decrease;  and  over  this  same  time  interval  the  drain 
current  decreases.  Similarly,  there  is  a  decreasing  gate  current  over  this  time  interval, 
suggesting  again  strong  coupling  between  the  gate  and  the  drain  regions. 

The  displacement  currents  through  the  source  and  drain  regions  are  displayed  in  figure 
4.5c.  The  displacements  currents  which  are  qualitatively  similar  to  that  discussed  in 
connection  with  sudden  changes  in  bias,  occur  when  the  voltage  changes  across  the  gate 
are  at  their  steepest  variation,  and  suggest  that  their  origin  is  in  the  oscillation  of  the 
charge  density  in  the  2DEG. 

In  the  calculations  associated  with  the  sudden  voltage  change  the  high  frequency 
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Figure  4.3  T ransient  Response  for  a  Change  in  Bias,  Vfo  =  1.5  to  4.0  and  VgS 
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Figure  4.4  T  ransient  Response  for  a  Change  in  Bias,  V^s  =  1.5  to  4.0  and  VgS  =  0.4 
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charge  density  variations  occured  in  the  2  DEG  within  the  source  and  drain  regions.  The 
variation  of  the  charge  under  the  gate  is  dominated  by  the  gate  contact.  For  the 
calculation  here,  we  display  the  density,  temperature,  magnitude  of  velocity,  and  potential 
variation  at  point  ’2’  of  figure  3.6.  The  steady  state  results,  which  incorporate  velocity 
overshoot,  show  an  additional  enhancement  of  velocity  overshoot,  along  with  substantial 
changes  in  the  temperature  of  the  carriers.  The  carrier  velocity,  while  not  periodic,  tends 
to  show  structure  at  the  end  of  the  10  ps  cycle,  indicating  that  is  closely  following  the 
changes  in  gate  bias.  The  electron  temperature  variation  does  not  display  a  full  two  cycles 
at  the  end  of  10  ps,  an  indication  that  the  relaxation  effects  require  a  longer  time  to  settle 
in.  After  some  point  in  time  all  of  the  relevant  quantities  such  as  energy,  density, 
temperature,  etc.,  will  be  periodic  with  a  period  of  5  ps.  At  this  point  the  temperature  will 
undergo  smaller  excursions.  While  the  oscillation  has  not  yet  settled  to  a  steady  state 
oscillation,  the  results  again  indicate  that  the  charge  variations  under  the  gate  are 
determined  by  the  time  dependence  of  the  gate  contact,  but  that  the  oscillation  in  the 
region  of  the  source  and  drain  contacts  are  significantly  modified  by  the  high  frequency 
charge  density  oscillations. 

The  principle  conclusion  of  the  transient  calculation  is  that  the  time  to  relaxation  of  the 
structure  is  dominated  by  the  longest  relaxation  time,  which  for  this  structure  and  material 
is  the  energy  relaxation  time;  and  that  this  relaxation  time  will  dominate  all  time 
dependent  oscillatory  behavior. 


TIME,  PSEC 

Figure  4.5a  Time  dependent  gate  voltage  on  the  PHEMT  for  a  fixed  drain  voltage  of  3.0  v. 
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Figure  4.5b  Total  current  out  of  all  three  contacts  following  the  time  dependent  imposition  of 
the  gate  voltage  of  figure  4.5  a 
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igure  4.5c  Displacement  current  out  of  the  source  and  drain  contacts  following  the  time 
tpendent  imposition  of  the  gate  voltage  of  figure  4.5a 
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Figure  4.5d  Density,  temperature,  velocity  magnitude,  and  potential  at  point  ’2’  following  the 
time  dependent  imposition  of  the  gate  voltage  of  figure  4.5a 


5.  Monte  Carlo  Simulations 


Introduction 

This  section  is  concerned  with  the  development  of  Monte  Carlo  procedures  for  examining 
transport  in  FETs.  The  basic  question  we  asked  is  why  perform  Monte  Carlo  calculations 
at  all;  and  more  specifically  why  perform  them  as  part  of  this  study.  There  are  several 
compelling  reasons:  First,  MC  procedures  offer  the  best  means  of  determining  the  exact 
form  of  the  distribution  function  for  semiclassical  structures.  The  moments  of  the 
Boltzmann  transport  equation  discussed  in  the  previous  sections  start  off  assuming  a 
specific  form  for  the  distribution  function.  Second,  when  transients  are  incorporated,  the 
exact  form  of  the  distribution  function  tends  to  reveal  transient  transport  effects  not 
apparent  from  the  moment  studies.  Third,  MC  procedures  offer  a  check  on  the 
assumptions  used  in  the  moment  formulation;  and  hence  provide  a  means  through  which 
the  scattering  parameters  and  dissipation  mechanisms  can  be  more  representatively 
handled  using  the  moment  methods.  This  third  reason  evolved  into  the  primary 
motivating  factor  of  incorporating  the  MC  procedures  in  the  study,  particularly  when  it 
was  recognized  that  the  transient  calculations  could  be  examined  more  efficiently  for 
complex  structures  using  moment  rather  than  MC  procedures.  All  calculations  were  for 
GaAs  using  standard  GaAs  parametric  values,  see,  e.g.,  Grubin  (1988). 

The  discussion  that  follows  describes  the  evolution  of  the  MC  development  and  its 
application  to  FETS.  All  calculations  performed  were  for  two  dimensional  structures; 
although  initially  the  two  dimensional  structure  was  a  two  terminal  diode  with  one 
dimensional  symmetry.  Performing  two  dimensional  MC  simulations  on  a  one 
dimensional  structure  permitted  us  to  expose  the  problems  associated  with  MC  methods 
without  introducing  complications  associated  with  two  dimensional  boundary  condition 
effects.  Parallel  MC  and  moment  calculations  were  performed  for  the  FET. 

Diode  Simulations 

As  discussed  above,  prior  to  applying  the  Monte  Carlo  simulation  procedure  to  a 
three-terminal  device,  such  as  an  FET,  the  procedure  was  applied  to  a  two-terminal, 
uniformly-doped  device.  The  device,  shown  in  Figure  5.1,  is  0.2  pm  wide  with  a  source 
drain  spacing  of  0.5  /im,  and  is  uniformly  doped  to  2  x  101 7 /cm*.  In  this  calculation  the 
component  ’z’  is  between  the  source  and  the  drain,  while  the  component  V  is  normal  to 
’z\  We  note  that  while  this  structure  would  yield  one-dimensional  solutions  using  a  drift  and 
diffusion  or  MBTE  approach,  the  Monte  Carlo  solution  may  exhibit  some  two 
dimensionality  in  the  solution  when  a  ” snapshot”  of  the  solution  is  taken  at  any  given  instant. 
However,  when  averaged  over  time  appropriately  the  one-dimensional  nature  of  the 
solution  is  revealed.  This  is  realized  as  we  examine  the  solutions  for  this  structure. 


-57- 


The  one  dimensional  calculations  are  summarized  by  current  voltage  characteristics  that 
are  nonlinear  and  reflect  the  nonuniform  distribution  of  carriers  in  the  structure.  Figure 
5.2  is  the  predicted  current-voltage  characteristic  for  the  structure.  (The  details  of  the 
current  calculation  are  discussed  in  the  FET  (next)  section.)  From  this  result  it  is 
apparent  that  the  low  field  mobility,  as  determined  from  the  slope  of  the  curve  near  the 
origin,  is  approximately  2,500  cm/volt-sec.  This  value  is  in  reasonable  agreement  with  that 
reported  by  others.  We  also  note  that  the  current  reaches  a  level  of  1160A/m  at  a  bias  of 
1.0  volt.  This  yields  a  mean  velocity  <v>  -  (j/e)/(2xl01 7/cm3)  =  1.8xl07 cm/sec  for  a 
mean  field  <F>  -  (1.0  v)/(0.5  n m)  =  20  kv/cm,  and  reflects  the  non-equilibrium  nature 
of  the  transport  within  the  device. 

Figure  5.3  shows  the  distribution  of  the  components  of  velocity  in  the  x  and  z  directions 
and  the  energy  level  of  carriers  in  the  gamma  valley  as  a  function  of  distance  (z)  between 
the  cathode  (z  =  0)  and  the  anode  (z  =  .5).  Consider  figure  5.3a  which  is  a 
representation  of  the  x  component  of  r  valley  velocity.  The  x  component  is  parallel  to  the 
contacts.  As  discussed  above  two  dimensional  calculations  were  performed  for  a  two 
terminal  device  with  one  dimensional  symmetry.  Thus  the  x  component  of  velocity  should 
fill  a  three  dimensional  array  with  side:  (vx,  x,  z)  as  shown  in  the  inset  of  figure  5.3a.  The 
distribution  of  points  in  figure  53a  represents  a  projection  of  all  points  in  the  MC 
simulation  onto  one  plane.  Each  point  in  this  figure  represents  one  particle,  and  it  is 
apparent  that  the  x  component  of  velocity  sustains  a  mean  value  of  zero.  While  the  mean 
x  component  of  velocity  is  zero,  we  observe  random  fluctuations  as  great  as  ±9.8  x 
107  cm/sec.  However,  the  majority  of  particles  are  in  the  range  of  ±2.8  x  107  cm/sec.  We 
also  note  a  lower  population  of  gamma  valley  particles  towards  the  anode  (z  =  0.5).  This 
is  due  to  transfer  to  the  upper  valley. 

Figure  5.3b  shows  a  projection  onto  one  plane  of  the  z  component  of  velocity.  Here  the 
mean  value  of  the  velocity  of  gamma  valley  carriers  is  near  5.6  x  107  cm/sec  in  the  middle 
of  the  device,  a  value  consistent  with  the  results  of  Grubin  and  Ferry  [  ]  for  a  field  near  20 
kv/cm.  Towards  the  anode  end  of  the  device  the  gamma  valley  carrier  velocity  is 
somewhat  higher,  but  fewer  carriers  remain  in  this  valley. 

Figure  5.3c  displays  a  projection  of  the  energy  distribution  for  the  gamma  valley  carriers. 
The  carriers  enter  in  equilibrium  with  the  lattice  and  gain  energy  under  the  influence  of 
the  applied  field,  remaining  within  the  gamma  valley  for  nearly  60%  of  the  structure.  The 
population  decreases  significantly  as  the  carriers  approach  the  anode. 
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Figure  5.1  Sketch  of  the  one  dimensional  diode  structure  used  in  the  Monte  Carlo  studies. 


Figure  5.2  Sketch  of  current  versus  voltage  for  the  structure  used  in  the  Monte  Carlo  studies. 
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Figure  5.3a  Projection  of  the 
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Figure  5.3c  Projection  of  the  distribution  of  T  energy  in  the  V  direction  for  a  bias  of  l.Ov. 


Figures  5.4  show  similar  plots  for  L  valley  carriers.  The  first  observation  in  all  of  the  plots 
shown  in  Figure  4  is  the  increasing  population  of  L  valley  carriers  as  the  anode  is 
approached.  The  x  component  of  velocity,  as  in  the  gamma  valley  result,  exhibits  a  mean 
of  zero.  The  z  component  of  velocity  exhibits  a  slight  positive  mean.  However,  this  value 
is  significantly  below  that  of  the  gamma  valley  carriers  due  to  the  greater  effective  mass  of 
L  valley  electrons.  Finally,  Figure  5.4c  shows  the  energy  distribution.  Note  the  large  L 
valley  energy  distribution  at  the  anode,  representing  the  consequences  of  electron 
transfer. 

As  was  indicated  in  the  discussion  previously,  in  which  we  are  solving  a  one  dimensional 
problem  within  a  two  dimensional  grid,  the  instantaneous  statistical  distribution  of 
electrons  and,  hence,  potential  within  the  device  displays  an  apparent  two  dimensional 
structure.  This  is  due  to  the  limited  number  of  particles  used  in  the  simulation,  and  as  the 
number  of  particles  increases,  the  deviation  from  the  one-dimensional  norm  should 
asymptotically  reach  zero.  To  show  this  two-dimensional  effect  the  net  charge  distribution 
along  lines  of  constant  x,  running  from  cathode  to  anode,  at  x  values  of  0.02,  .06,  .10,  .14, 
and  .18  ^m  are  presented  in  Figure  5.5. 
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Figure  5.4a  Projection  of  the  distribution  of  L  velocity  in  the  V  direction  for  a  bias  of  l.Ov. 
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Figure  5.4b  Projection  of  the  distribution  of  L  velocity  in  the  ’z’ direction  for  a  bias  of  l.Ov. 
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Figure  5.4c  Projection  of  the  distribution  of  L  energy  in  the  V  direction  for  a  bias  of  l.Ov. 

The  net  charge  distribution  at  point  (z)  is  -e(n-Nj)),  and  if  the  density  were  constant  and 
equal  to  background  the  net  charge  distribution  would  be  centered  about  zero.  For  the 
one  dimensional  diode  problem,  the  initial  condition  of  net  charge  neutrality  was  imposed. 
For  this  solution,  which  yields  a  nonuniform  potential  distribution  that  increases  from 
cathode  to  anode,  there  is  a  subsequent  accumulation  of  carriers  in  the  vicinity  of  the 
anode  and  a  depletion  of  carriers  in  the  vicinity  of  the  cathode. 

As  seen  in  figure  5.5  the  net  charge  distributions  are  negative  in  the  vicinity  of  the  cathode 
(indicating  charge  depletion)  and  positive  in  the  vicinity  of  the  anode  (indicating  charge 
accumulation).  The  charge  distributions  in  figure  5.5  are  obtained  from  a  cloud-in-cell 
smoothing  of  the  particles  for  use  in  Poisson’s  equation.  The  particular  results  shown  are 
for  a  1  volt  bias  across  the  device.  As  can  be  seen,  there  is  considerable  scatter  in  the 
results,  but  all  plots  show  some  degree  of  depletion  at  the  cathode  and  accumulation  at  the 
anode.  We  note  that  at  the  cathode,  electrons  are  injected  to  maintain  charge  neutrality. 

The  potential  variation  along  the  same  planes  in  the  device  are  shown  in  Figure  5.6.  The 
variation  across  the  device,  while  still  apparent,  is  less  obvious  due  to  the  smoothing 
nature  of  solutions  to  Poisson’s  equation.  The  curvature  of  the  resulting  distribution  of 
potential  reflects  the  depletion  and  accumulation  of  carriers  within  the  structure.  As  is 
apparent,  the  solution  is  not  a  uniform  field  solution. 
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Figure  5.5  Normalized  net  charge  distribution  along  the  z  direction  at  five  different  values  of 
x :  0.02,  0.06,  0.10,  0.14,  0.18  microns. 
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Figure  5.6  Potential  distribution  along  the  z  direction  at  five  different  values  of  x :  0.02,  0.06, 
0.10,  0.14,  0.18  microns. 
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GaAs  FET  Svrudations 


Test  Structure,  VgS  =  0.1  v,  V^s  =  0.5  v 

The  first  structure  considered,  shown  in  Figure  5.7,  was  a  small  modification  of  the  diode 
structure,  and  enabled  us  to  develop  procedures  for  studying  FETs  using  MC  methods. 
The  device  length  is  0.6  tim,  with  a  centrally  placed  0.2  ^m  gate.  The  channel  depth  was 
taken  as  0.2  um  and  uniformly  doped  to  1  x  1017/cm3.  This  device  constitutes  a 
normally-on  device.  The  Schottky  gate  depletion  layer  thickness  was  estimated  to  be  0.1 
um.  A  simulation  was  run  for  a  drain  bias,  V^g  =  0.5  volts  and  a  gate  bias,  vgs  '  'U 
volts.  Four  hundred  time  steps  (Poisson  updates)  were  used. 

Calculation  of  Current 

For  these  calculations,  which  were  taken  to  steady  state  it  was  assumed  that  no  particle 
current  passed  through  the  gate  contact;  reflecting  boundary  conditions  at  the  gate  were 
imposed.  The  source  and  drain  currents  were  calculated  by  determining  the  numbers  of 
particles  passing  through  the  source  and  drain  regions  as  a  function  of  time,  figure  5.8.  In 
steady  state  the  curves  are  linear  with  the  slope  (drawn  in  each  of  these  figures)  of  these 
curves  yielding  the  current  at  the  contacts.  Note  the  nonlinear  time  dependence  of  the 
particle  history  at  the  early  time  stages.  For  this  calculation  the  gate  particle  current  is 
zero  and  the  source  and  drain  currents,  agree  to  within  numerical  accuracy;  the  source 
current  is  5.41  x  102A/m,  the  drain  current  is  5.49  x  102A/m.  These  currents  are 
approximately  a  factor  of  two  lower  than  the  current  obtained  with  the  simple  diode 
structure  discussed  earlier;  and  may  be  accounted  for,  in  part,  due  to  differences  in  the 
background  charge  density,  and  under  the  gate  by  much  higher  values  of  velocity, 
indicating  the  prominence  of  velocity  overshoot,  as  discussed  below. 
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Figure  5. 7  Sketch  of  the  GaAs  FET  structure  used  in  the  Monte  Carlo  studies. 
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Figure  5.8aTime  history  of  charge  across  the  source  contact  for  Vq  =  !  v,  Vjy-O.Sv 


Figure  5.8b  Time  history  of  charge  across  the  drain  contact  for  Vq- 0.1  v,  V[)  =  0.5  v 
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Figure  5.9  Distibution  of  charge  at  the  termination  of  the  calculation  for  Vq  =  0.1  v, 
Vj)  =  0.5  v 

The  distribution  of  all  particles  in  the  device  at  the  termination  of  the  calculation  is  shown 
in  Figure  5.9.  Here  we  observe  the  expected  result  of  a  uniform  distribution  of  electrons 
under  the  source  and  drain  with  the  depletion  under  the  gate  clearly  evident. 

Figure  5.10  shows  the  projections  of  the  distributions  of:  (a)  r  particle  velocity  along  the  z 
direction  parallel  to  the  channel  (figure  5.10a);  (b)  r  particle  velocity  along  the  x 
direction  normal  to  the  channel  (figure  5.10b);  and  (c)  T  energy  (figure  5.10c).  The 
results  are  only  apparently  similar  to  the  diode  calculation,  because  the  current  paths 
through  the  device  suggest  that  the  higher  carrier  velocities  will  not  be  uniformly 
distributed  along  the  channel.  But  we  note  that  the  calculations  show  a  velocity 
distribution  in  the  vicinity  of  the  source  contact  similar  to  that  associated  with  the  diode 
calculation  (figure  5.3b);  this  is  expected  as  the  current  density  is  a  factor  of  two  smaller, 
and  the  background  density  is  a  factor  of  two  smaller  than  the  diode  calculation.  Under 
the  gate  region  the  velocity  distribution  displays  larger  values  than  the  diode  calculation, 
as  expected  on  the  basis  of  current  continuity. 

Figuure  5.10b  displays  the  expected  result  that  the  mean  component  of  the  x  component  of 
r  velocity  is  near  zero.  The  energy  variation  (figure  5.10c)  shows  the  greatest  increase  at 
the  drain  side  of  the  gate,  and  while  there  are  still  high  energy  carriers  at  the  drain,  there 
are  also  a  greater  number  of  lower  energy  carriers  there  as  electrons  relax  in  the  lower 
field  at  the  drain.  The  results  for  the  L  valley  electrons  are  shown  in  Figure  5.11.  As  the 
figures  show,  by  the  low  number  of  particles,  very  few  electrons  transfer  to  the  L  valley 
and  their  velocity  mean  is  near  zero  in  both  the  x  and  z  directions. 
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Figure  5.10c  Projection  of  the  distribution  of  r  energy  in  the  ’z’ direction. 
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Figure  5.11a  Projection  of  the  distribution  of  L  velocity  in  the  V  direction. 
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Figure  5.11b  Projection  of  the  distribution  of  L  velocity  in  the  ’x’  direction. 


Figure  5.11c  Projection  of  the  distribution  of  L  energy  in  the  ’z’ direction. 
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FET  Simulations,  Monte  Carlo  and  Moment  Equation  Comparisons 

Full  MC  FET  simulations  were  performed  for  the  structure  shown  in  figure  5.7,  but  with 
the  channel  width  reduced  to  0.1  pm.  These  calculations  were  then  compared  to  parallel 
calculations  performed  using  the  moments  of  the  Boltzmann  transport  equation. 

A  complete  set  of  current-voltage  curves  was  generated  for  this  device.  These 
characteristics  are  shown  in  Figure  5.12.  (An  additional  set  of  calculations  were  performed 
for  a  similar  device  with  a  0.1  pm  undoped  substrate  at  VgS  =  0.3  volts.  These  calculations 
yield  a  current  level  approximately  four  times  higher  over  a  large  bias  range  that  for  the 
basic  0.1  pm  channel  height  structure  at  VgS  =  0.3  v.)  The  results  for  the  device  without  a 
substrate  are  in  good  agreement  with  those  reported  for  a  similar  device  others.  For  this 
device,  without  the  substrate,  the  drain  current  at  VgS  =  0.1  and  =  0.5  volts  is  about 
16A/m  or  about  35  times  smaller  than  the  0.2  ^im  channel  depth  device  discussed 
previously;  this  device  is  ’normally-off . 


Figure  5.12  Sketch  of  current  versus  voltage  for  the  0.1  micron  channel  height  FET. 
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Detailed  Results  for  V^s  =  0.5  v,  VgS=  0.1  v. 


The  time  history  of  the  charge  crossing  the  source  and  drain  contacts  for  VgS=  0.1, 
Vds  =  0.5  is  shown  in  Figure  5.13.  While  the  scales  for  the  source  and  drain  charge  are 
different,  the  calculations  indicate  that  it  takes  a  longer  time  for  the  source  contact  to 
reach  steady  state.  At  the  end  of  10  psec,  the  calculation  was  terminated  with  the  resulting 
slopes  at  the  source  and  drain  yielding  current  levels  at  the  source  of  15.4  A/m,  and  at  the 
drain  of  18.5  A/m. 

The  total  electron  distribution  in  the  device  at  V^s=  0.5  v  and  VgS  =  0.1  v  is  shown  in 
Figure  5.14.  Here  we  observe  the  complete  depletion  of  carriers  under  the  gate,  pinching 
off  the  channel  and  an  extension  of  the  depletion  region  downstream  of  the  gate  towards 
the  drain. 

Figure  5.15  shows  the  projection  of  the  distribution  of  the  velocity  components  and  energy 
for  the  gamma  valley  carriers  at  this  bias.  Both  the  z  component  of  velocity  parallel  to  the 
channel  (figure  5.15a),  and  the  x  component  of  velocity  normal  to  the  channel  (figure 
5.15b),  are  very  small.  From  the  current  it  was  determined  that  the  mean  velocity  of 
carriers  crossing  the  source  and  drain  contacts  is  of  the  order  of  10®  cm/sec.  We  note  that 
the  few  gamma  valley  carriers  under  the  gate  have  a  much  higher  velocity  ranging  from  3  x 
107  to  6  x  107  cm/sec.  Recall  that  these  projection  plots  are  instantaneous  distributions  at 
the  end  of  the  computation  and  do  not  reflect  the  variation  which  may  be  expected  from 
time  increment  to  time  increment.  Due  to  the  low  velocity  at  the  source  and  drain  and  low 
current,  the  gamma  valley  carrier  energy  levels  are  very  low,  as  shown  in  Figure  5.15c. 
The  results  for  the  L  valley  carriers  are  not  shown  for  this  case,  since  very  few  carriers 
transfer  from  the  gamma  valley. 

Detailed  Results  for  V^s  —  0.5  v,  VgS  =  0.5  v. 

When  the  gate  bias  level  is  increased  to  VgS  =  0.5  volts,  while  holding  at  0.5  volts,  the 
electron  distribution  is  as  shown  in  Figure  5.16.  There  are  significantly  more  particles  in 
the  channel  under  the  gate  at  this  bias  level,  indicating  the  reduction  in  the  depletion  layer 
thickness.  The  current  level  is  a  higher  160A/m. 

Figure  5.17  displays  the  projections  of  the  z  component  of  gamma  velocity  (figure  5.17a), 
the  x  component  of  gamma  velocity  (figure  5.17b),  and  the  gamma  energy  (figure  5.17c). 
When  compared  to  Figure  5.15,  we  observe  substantially  more  high  velocity  carriers 
throughout  the  structure,  and  a  correspondingly  broader  distribution  of  higher  energy 
carriers  in  figure  5.17.  There  are  also  many  more  carriers  in  the  device.  This  occurs 
because  of  the  greater  depletion  in  the  VgS  =  0.1  volt  case. 
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Figure  5.14  Distibution  of  charge  at  the  termination  of  the  calculation  for  Vq—0.1  v, 
Vd  =  0.5v 


Figure  5.15a  Projection  of  the  distribution  of  T  velocity  in  the  ’z’  direction  for  Vq  =  0.1, 
Vd  =  0.5v  _ 
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Figure  5.15b  Projection  of  the  distribution  of  r  velocity  in  the  V  direction. 


Figure  5.15c  Projection  of  the  distribution  of  r  energy  in  the  ’z’  direction. 


Several  additional  generic  comments  are  in  order:  (i)  In  the  MC  simulation,  3,000 
particles  represent  a  charge  neutral  situation.  As  the  depletion  layer  evolves,  the  total 
number  of  particles  in  the  device  drops  below  3000  to  a  number  dependent  on  the  gate 
bias  and  depletion  layer  structure,  (ii)  Comparing  figures  5.15b  and  5.17b  we  recall  that 
the  mean  z  velocity  component  at  the  source  and  drain  in  Figure  5.15b,  where  the  current 
level  was  16A/m,  was  10®  cm/sec.  In  Figure  5.17b,  where  the  current  level  is  160A/m,  the 
mean  velocity  at  the  source  and  drain  is  107  cm/sec.  This  result  is  difficult  to  obtain  from 
the  velocity  projections.  It  is  however  apparent  that  very  high  particle  velocities, 
approaching  9  x  107  cm/sec  and  above,  are  present  in  the  channel,  (iii)  The  energy 
distribution,  shown  in  Figure  5.17c,  shows  energy  levels  in  excess  of  0.3  eV,  and  electron 
transfer  would  be  expected. 

The  L  valley  carrier  velocity  and  energy  distributions  are  shown  in  Figures  5.18a-c.  The 
figures  indicate  transfer  occurs  at  the  drain  end  of  the  device,  somewhat  past  the  end  of 
the  gate.  The  carrier  velocities  remain  low,  as  does  the  L  valley  energy,  due  to  the  higher 
effective  mass.  The  point  to  note  is  that  while  the  numbers  of  particles  transferred  to  the 
L  valley  is  substantial,  current  is  still  dominatedy  the  fewer  gamma  valley  carriers  with 
their  higher  velocities. 
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Figure  5.16  Distribution  of  charge  at  the  termination  of  the  calculation  for  Vq  =  0.5  v, 
VD-0.5  v 
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Figure  5.17a  Projection  of  the  distribution  of  T  velocity  in  the  'z'  direction  for  Vq  -  0.5, 
Vd  =  0.5v 


Figure  5. 1 7b  Projection  of  the  distribution  of  T  velocity  in  the  ’x’  direction. 


Figure  5.17c  Projection  of  the  distribution  of  V  energy  in  the  ’z’  direction. 


Figure  5.18a  Projection  of  the  distribution  of  L  velocity  in  the  ’z’  direction  for  Vq  =  0.5, 
V D~  0.5  v. 


Figure  5.18b  Projection  of  the  distribution  oj  !.  velocity  in  the  V  direction. 
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Figure  5.18c  Projection  of  the  distribution  of  L  energy  in  the  ’z’  direction. 
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Comparison  to  MBTE  Calculations 


The  MC  calculations  were  compared  to  calculations  using  the  moments  of  the  Boltzmann 
transport  equation  algorithm.  This  was  the  algorithm  implemented  for  the  study  of  the 
pseudomorphic  HEMT.  For  this  comparison  a  complete  set  of  current  voltage 
characteristics  was  obtained.  As  discussed  in  the  earlier  chapters  the  calculations 
performed  with  the  moment  equations  involve  the  incorporation  of  scattering  integrals 
whose  value  depends  upon  assume  scattering  constants  such  as  the  acoustic  deformation 
coupling  constant,  the  Frohlich  field,  etc.  Very  often  in  performing  these  calculations  the 
values  chosen  for  the  parameters  are  adjusted  until  a  uniform  field  dependent  velocity  is 
computed  and  compared  to  a  similar  calculation  using  MC  procedures.  Adjustments  in 
the  MBTE  field  dependent  velocity  parameters  are  then  made  until  agreement  with  the 
MC  result-  is  obtained. 

Two  comparative  calculations  were  performed  here.  In  the  MB!  E  calculation  only  two 
velleys  were  considered.  In  one  case  the  r  -L  separtion  was  0.33  ev,  which  was  greater  than 
the  0.30  value  used  in  the  MC  study.  In  the  second,  the  r-L  separation  was  0.284  ev. 
Additionally,  in  the  larger  energy  separation  calculation  the  static  dielectric  constant  was 
taken  as  12.53.  In  the  lower  intervalley  calculation  the  static  dielectric  constant  was  set  to 
13.81.  In  all  three  calculations  the  current  voltage  characteristics  for  VgS  =  0.1  v,  were 
similar  as  shown  below. 

Before  discussing  the  current  voltage  characteristics  it  is  useful  to  compare  the  detailed 
results  for  the  MC  and  MBTE  calculations  at  a  specific  bias  point:  =  1.0  v,  VgS  =  0.3 

v.  For  example  figure  5.19  is  a  comparison  of  the  density  distribution  obtained  from  both 
the  MC  and  the  MBTE  calculations.  The  top  frame  displays  contours  of  constant  charge 
density.  The  half-ellipses  represent  the  shape  of  the  depletion  region.  It  appears  to  be 
qualitatively  similar  to  the  shape  of  the  depletion  region  obtained  from  the  MC 
calculation. 

Figure  5.20  is  a  projection  of  the  potential  surface  for  this  calculation  at  =  1.0  v  and 
VgS  =  0.3  v.  Note  that  for  the  MBTE  calculation  and  on  the  channel  side,  the  potential 
increases  monotonically  from  the  source  to  a  value  one  volt  higher  at  the  drain  contact. 
The  dip  in  potential  on  the  gate  side  represents  the  effect  of  the  gate  boundary  condition. 
For  the  MC  calculations,  apart  from  some  irregularities  in  the  surface,  due  to  the 
statistical  nature  of  the  MC  particle  distribution,  the  potential  surfaces  are  very  similar. 

The  predicted  current-voltage  characteristics  for  this  structure  as  obtained  from  the  MC 
and  MBTE  are  compared  in  figures  5.21  and  5.22.  In  figure  5.21,  the  larger  intervalley 
energy  separation  was  assumed,  in  figure  5.22,  the  smaller  value  was  used.  The  curves 
labeled  "A”  were  obtained  for  a  finite  but  small  value  of  the  thermal  conductivity  (see 
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section  2)  As  a  result,  the  heating  of  the  r  valley  carriers  is  a  maximum  and  remains 
highly  localized  due  to  minimal  heat  conduction.  This  gives  the  highest  level  of  intervalley 
transfer  and  results  in  the  observed  negative  conductance  at  higher  gate  voltages,  which  is 
in  disagreement  with  the  MC  result.  The  cause  of  the  discrepancy  lies  in  either  the  value 
of  the  thermal  conductivity  chosen,  or  the  scattering  rates,  or  both.  That  is,  the  predicted 
r  valley  temperature  is  too  high,  or  at  the  predicted  temperatures  the  intervalley  scattering 
is  too  large.  We  emphasize  that  under  uniform  field  conditions  both  the  MC  and  the 
MBTE  approaches  yield  similar  field  dependent  velocity  curves.  Thus  the  result  is  further 
complicated  by  non-uniform  field  effects.  If  the  problem  lies  with  the  predicted  levels  of 
the  r  valley  temperature,  then  the  thermal  conductivity  and  the  momentum  and  energy 
relaxation  times  may  be  suspect.  Assuming  the  problem  is  limited  to  thermal  conductivity, 
a  limit  may  be  explored  by  setting  the  thermal  conductivity  to  infinity.  Under  such 
conditions  the  electron  temperature  will  remain  constant  everywhere  within  the  structure  at 
the  lattice  temperature  (electrons  enter  the  device  at  equilibrium  with  the  lattice  and  the 
electron  temperature  reflects  this.)  Such  a  condition  results  in  the  current  voltage 
characteristics  shown  as  "B"  in  figure  5.21.  From  this  result  we  observe  that  the 
manipulation  of  the  thermal  conductivity  brackets  the  MC  result  and  the  numerical  results 
could  be  matched  if  desired.  However,  this  approach  is  not  satisfactory. 

A  second  investigation  was  performed  using  the  initial  finite  value  of  thermal  conductivity 
but  with  the  smaller  value  of  r-L  separation.  The  field  dependent  velocity  is  nearly  the 
same  for  both.  The  MBTE  calculation  performed  using  scattering  rates  based  on  these 
parameters  yield  the  current  voltage  characteristics  presented  in  figure  5.22.  Here  we  see 
excellent  agreement  between  the  MC  and  MBTE  results. 

Without  editorializing  about  better  or  worse  agreement  between  different  sets  of 
numerical  procedures  and  different  sets  of  calculations,  the  results  tend  to  isolate 
differences  in  the  physical  assumptions  and  consequences  of  the  two  approaches.  We  note 
that  in  the  MC  calculation  carriers  with  an  energy  sufficient  for  intervalley  transfer 
undergo  such  transfer.  In  the  MBTE  where  we  are  dealing  with  mean  energies,  at  every 
value  of  energy  there  are  always  carriers  in  the  upper  valley.  Thus  the  transition  is  not 
sudden.  If  the  r-L  separation  is  larger  than  that  used  for  the  MC,  then  there  is  a  delay  in 
the  field  value  at  which  intervalley  transfer  occurs.  This  delay  in  field  value  does  not 
imply  that  negative  conductance  in  the  IV  characteristics  will  necessarily  occur,  for  the 
rate  of  intervalley  transfer  is  controlled  by  a  number  of  parameters  including  the  coupling 
constants.  The  best  that  can  be  said  in  comparing  the  MC  and  MBTE  calculations  is  that: 
(a)  they  both  provide  qualitatively  similar  physics,  and  (b)  adjustments  of  parameters  can 
lead  to  close  quantitative  numerical  agreement.  The  choice  of  using  either  MC  or  MBTE 
then  rests  with  other  issues,  such  as  ease  of  the  compuxational  procedure,  etc.! 
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Figure  5.21  Comparison  of  I-V  Characteristics  as  determined  from  the  Monte  Carlo,  and 
MBTE  solutions  with  Ap.jr^  =  0.36  ev. 


Figure  5.22  Comparison  of  I-V  Characteristics  as  determined  from  the  Monte  Carlo,  and 
MBT E  solutions  with  A —  0.284  ev. 
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6.  One  Dimensional  Transient  Simulations:  Search  for  Overshoot 


The  search  for  velocity  overshoot  has  occupied  the  attention  of  many  workers  over  the  last 
decade.  The  difficulty  lies  in  measuring  velocity  overshoot,  since  it  is  always  present. 
Furthermore  all  electrical  measurements  involve  the  time  dependent  passage  of  charge, 
which  in  turn  is  sensitively  affected  by  local  strong  and  weak  nonuniformities  in  the 
distribution  of  carriers.  Additionally,  as  the  time  scales  are  reduced,  a  large  number  of 
calculations,  including  those  discussed  in  section  4,  indicate  that  the  transients  associated 
with  displacement  currents  occur  on  time  scales  commonly  associated  with  velocity 
overshoot  transients.  Nonetheless,  both  transient  and  steady  state  experiments  have  been 
deviced  to  identify  velocity  overshoot.  For  example:  (i)  In  the  early  70’s  it  was 
demonstrated  that  the  transient  operation  of  GaAs  would  exhibit  microwave  heating  and 
electrical  behavior  that  would  be  directly  attributable  to  velocity  overshoot,  while  (ii) 
Laval  et  al  (1980)  attempted  to  show  that  the  current  voltage  curves  would  show  an 
increase  in  current  saturation  as  the  active  region  of  the  structure  was  decreased  and  that 
this  increase  in  current  saturation  could  be  attributed  to  velocity  overshoot.  None  of  these 
approaches  is  regarded  by  the  present  authors  as  representative  of  velocity  overshoot. 
Rather  we  should  not  be  looking  at  velocity  overshoot  but  at  the  more  general  problem  of 
nonequilibrium,  since  during  the  course  of  a  transient  there  are  essentially  four  time 
scales  to  deal  with:  (i)  carrier  relaxation,  (ii)  momentum  relaxation,  (iii)  energy  relaxation, 
and  (iv)  displacement  current  (capacitive),  time  scales. 

The  simplest  type  of  experiment  to  exhibit  velocity  overshoot  is  a  two  terminal  experiment 
that  isolates  the  different  time  scales  associated  with  energy,  momentum  and  capacitive 
contributions.  In  section  4,  involving  the  transient  behavior  of  the  pseudomorphic 
HEMT,  it  was  determined  that  the  time  to  relaxation  was  determined  by  the  longest  time 
scale,  which  was  the  relaxation  of  the  electron  temperature.  The  two  terminal  transient 
simulations  discussed  below  are  intended  to  highlight  the  differences  between  the 
temperature,  momentum  and  capacitive  time  scales-and  to  identify  the  differences  as 
representative  of  nonequilibrium  phenomena. 

As  discussed  by  the  present  workers,  the  issue  of  velocity  overshoot  is  nonlinear  in  the 
sense  that  at  low  values  of  bias  there  is  very  little  electron  transfer,  at  very  high  values  of 
bias  there  is  a  great  deal  of  electron  transfer.  If  the  signals  imposed  on  a  device  do  not 
result  in  the  transfer  of  carriers  from  the  lower  to  upper  valleys  or  from  the  upper  to  lower 
valleys,  then  any  nonequilibrium  contributions  arise  from  differences  between  energy  and 
momentum  relaxation  time  without  the  competing  effects  of  intervalley  transfer.  At 
intermediate  values  of  bias  intervalley  transfer  occurs  and  complicates  the  issue.  To 
determine  the  means  of  resolving  some  of  these  competing  effects  we  consider  the 
following  numerical  experiments.  The  transient  calculations  are  obtained  from  a  solution 
to  the  moments  of  the  Boltzmann  transport  equation  for  the  structure  of  figure  6.1. 
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Figure  6.1  Structure  used  in  the  transient  simulations. 


We  consider  the  NIN  structure  identified  in  figure  6.1.  The  end  regions  are  doped  to 
10 1 7  /cm3  with  the  interior  2.0  micron  region  being  undoped.  The  structure  is  subject  to  a 
DC  bias  of  2.0  v,  and  reaches  steady  state  with  a  current  level  of  4.2A/cm2 . 

Superimposed  upon  the  DC  result  is  a  time  dependent  AC  voltage.  The  time  dependent 
voltage  is  periodic  exhibiting,  repetitive  excitation  and  relaxation.  The  excitation  portion 
generally  results  in  an  increase  in  current  beyond  the  dc  value.  The  relaxation  phase 
results  in  a  decrease  in  current  toward  the  steady  state  value.  As  a  baseline  calculation 
consider  figure  6.2. 

Figure  6.2a  displays  the  applied  voltage,  and  the  resulting  total  current  response  for  the 
diode.  The  voltage  train  consists  of  a  step  change  in  voltage  from  a  DC  value  of  2.0  v,  to 
4.0  v,  followed  by  an  exponential  decay  to  2.0  v.  The  period  of  the  oscillation  is  5  ps  and  is 
of  the  order  of  the  relaxation  time  for  the  electron  temperature  as  discussed  in  section  5. 
The  numbers  listed  in  the  voltage  pulse  of  figure  6.2a  are  keyed  to  the  spatial  dependent 
distributions  in  the  remaining  parts  of  figure  6.2.  We  point  out  that  these  calculations 
were  for  the  first  ten  picoseconds  of  the  transient.  In  steady  state  all  of  the  variables 
computed  will  reach  steady  state  with  the  period  of  five  picoseconds.  As  will  be  apparent 
from  t’ie  discussion,  not  all  variables  reach  steady  state  during  the  first  ten  picoseconds. 
The  consequences  of  this  will  be  discussed. 
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The  current  oscillation  displayed  in  figure  6.2a  shows  characteristics  of  the  displacement 
current  contributions  identified  in  section  4.  Anticipating  later  results  (see  figure  6.5)  the 
high  frequency  oscillations  are  dependent  upon  the  doping  level  of  the  N  +  regions  and 
decay  more  rapidly  with  increasing  doping  density.  Recall  that  the  dielectric  relaxation 
time  decreases  with  increased  doping  level.  A  comparison  of  the  two  figures  indicates  the 
more  rapid  oscillation  and  decay  for  the  higher  doping  density  calculation.  Note  also  that 
the  steady  state  current  for  both  cycles  of  the  figure  o.2a  calculations  are  the  same.  For 
this  transient  there  is  little  that  can  be  extracted  concerning  nonequilibrium  effects. 
Nevertheless  the  details  of  this  calculation  are  worthwhile  considering. 

Figure  6.2b  is  a  plot  of  the  self  consistent  voltage  within  the  structure  during  one  cycle  of 
the  oscillation.  The  numbers  in  each  frame  correspond  to  the  numbers  on  the  oscillatory 
profile  of  figure  6.2a.  The  general  feature  of  the  potential  profile  is  that  there  is  a  slight 
retarding  potential  in  the  vicinity  of  the  cathode  N  + 1  interface  with  a  monotonic  increase 
in  potential  toward  the  anode  boundary  where  the  potentials  is  flat  at  the  anode  N + 
region.  The  distribution  of  charge  accompanying  this  potential  change  is  shown  in  figure 
6.2c,  with  breakdown  between  the  r  and  L  valley  contributions  represented  by  their 
magnitudes.  For  this  oscillation  most  of  the  carriers  are  in  the  r  valley,  with  most  of  the 
transfer  occuring  for  this  range  of  bias  near  the  downstream  N  +  region.  Nonequilibrium 
electron  transfer  occurs  near  the  downstream  interface.  The  gamma  valley  carrier  velocity 
corresponding  to  the  voltage  changes  is  displayed  in  figure  6.2d.  Within  the  undoped 
region  the  carriers  are  reaching  speeds  near  4xl07  cm/sec,  (without  any  apparent  way  to 
detect  them,  since  an  independent  means  of  ascertaining  the  exact  density  distribution  is 
not  possible.)  The  velocity  distribution  shows  complete  relaxation  during  the  latter  half  of 
the  cycle,  a  feature  consistent  with  the  short  momentum  relaxation  time.  The  gamma 
valley  electron  temperature,  shown  in  figure  6.2e  requires  a  little  more  time  to  relax,  but 
for  all  practical  purposes  the  structure,  has  relaxed  to  its  steady  state  value  within  the  time 
frame  of  five  picoseconds. 

The  calculations  of  figure  6.2  are  such  that  the  structure  spends  little  time  in  the  high 
voltage  range.  There  is  a  sudden  increase  in  voltage  followed  by  an  exponential  drop  to  a 
steady  state.  Thus  the  detailed  transient  behavior  or  transfer  from  the  lower  to  the  upper 
valleys  and  then  return  is  not  exposed.  To  expose  this  detail  the  applied  voltage  must 
permit  some  fraction  of  the  electrons  to  undergo  transfer  to  and  from  a  specific  value.  We 
point  out  that  when  dealing  with  the  steady  state  equilibrium  velocity  field  curves  of  Phase  / 
the  issue  of  the  transient  transfer,  which  is  at  the  core  of  velocity  overshoot,  did  not  enter  the 
discussion.  Figure  6.3  displays  results  in  which  the  voltage  is  exponentially  increased  from 
two  volts  to  near  four  volts  and  exponentially  decreased  to  two  volts.  The  period  of 
oscillation  is  again  five  picoseconds.  The  keying  on  the  voltage  profiles  has  the  same 
significance  as  the  keying  of  figure  6.2. 


-87- 


CURRENT,  J/JREF  VOLTAGE 


0.00  2.00  4.00  6.00  8.00  10.00 

TIME,  PSEC 


Figure  6.2a  Voltage  and  current  profiles  for  the  nominal  NIN  structure,  with  a  step  change 
in  voltage  followed  by  an  exponential  decay.  The  dc  values  of  current  and  voltage  for  this 
calculation  are  V(DC)  =2.0v,  J  (DC)  =4.22x10*  A/cm2 . 
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Figure  6.2  b  Time  sequence  of  potential  proiles  corresponding  to  the  time  dependent 
calculation  of  figure  6.2a. 


Figure  6.2c  Time  sequence  of  density  proiles  corresponding  to  the  time  dependent 
calculation  of  figure  6.2a. 
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Figure  6.2d  Tittle  sequence  of  gamma  valley  velocity  profiles  corresponding  to  the  time 
dependent  calculation  of  figure  6.2a. 
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Figure  6.2e  Time  sequence  of  gamma  valley  temperature  profiles  corresponding  to  the  time 
dependent  calculation  of  figure  6.2a. 


There  are  several  important  features  to  note  about  the  figure  6.3  results  when  compared 
to  the  figure  6.2  results.  First ,  the  displacement  current  contributes  to  the  diode  current 
twice  during  the  oscillation,  once  as  the  voltage  increases  and  a  second  time  as  the  voltage 
decreases.  In  other  words  a  significant  change  in  the  structure  of  the  applied  voltage 
profile  introduces  the  displacement  current  effects.  Second,  the  current  level  to  which  the 
system  is  relaxing  to  after  the  displacement  current  effects  damp  out  is  different  and 
higher  for  the  figure  6.3  calculations  than  for  the  figure  6.2  calculations  (the  current  is 
negative,  and  higher  current  is  more  negative  current).  The  implication  is  that  momentum 
relaxation  has  occurred,  but  energy  relaxation  has  not  occurred.  Third,  on  the  first  cycle, 
the  steady  state  current  during  the  first  half  cycle  is  greater  than  the  steady  state  current 
during  the  second  half,  indicating  that  for  this  value  of  bias  transfer  from  lower  to  upper 
valleys  occurs  on  a  different  time  scale  than  transfer  from  the  upper  to  the  lower  valleys. 
During  the  second  part  of  the  cycle  this  difference  is  negligible,  as  indicated  by  the  inset, 
but  as  later  calculations  show  this  difference  is  prominent  at  other  voltage  levels. 

The  time  sequence  of  the  voltage  profiles  are  shown  in  figure  63b,  where  we  see  that  the 
voltage  drop  across  the  device  spends  a  longer  period  of  time  at  the  higher  values  than  the 
corresponding  profiles  of  figure  6.2b.  The  density  profiles  in  figure  6.3c  show  more 
electron  transfer  during  the  first  part  of  the  cycle  than  for  the  corresponding  profiles  of 
figure  6.2c.  The  gamma  valley  velocity  profiles  are  shown  in  figure  6.3d.  Note  that  for 
every  frame  the  velocity  is  higher  than  that  of  the  corresponding  figure  6.2d  profiles.  It  is 
this  higher  velocity  that  is  responsible  for  the  larger  value  of  steady  state  current.  The 
temperature  profiles  are  displayed  in  figure  63e.  We  note  that  everywhere  within  the 
structure  the  electron  temperature  is  higher  than  for  the  figure  6.2e  calculation.  At  the 
downstream  N  +  region  the  electron  temperature  is  highest  and  is  responsible  for  most  of 
the  electron  transfer.  It  is  interesting  to  point  out  that  at  the  end  of  each  cycle  for  the 
figure  6.3  calculation,  although  the  electron  temperature  is  higher  than  for  the  figure  6.2 
calculations,  more  of  the  diode  is  within  the  high  temperature  region  in  the  latter 
calculations.  > 

The  above  results  indicate  that  when  a  signal  with  a  period  of  the  order  of  the  energy 
relaxation  time  is  imposed  upon  the  structure,  the  detailed  intra-signal  shape  can 
determine  when  strong  nonequilibrium  effects  will  occur.  Differences  in  the  shape  of  the 
pulse  will  manifest  themselves  in  significant  differences  in  the  current  level  which  are  a 
direct  measure  of  the  differences  in  the  energy  and  momentum  relaxation  time. 
Essentially  the  calculations  of  figure  6.3  teaches  that  if  there  are  time  constants  associated 
with  the  intra-pulse  shape  and  these  time  constants  approach  the  momentum  relaxation 
time  then  the  current  levels  will  reflect  these  time  constants  through  increased  values.  The 
degree  to  which  these  relaxation  effects  appear  will  depend  when  the  intervalley  exchange 
is  most  prominent. 
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Figure  6.3a  Voltage  and  current  profiles  for  the  nominal  NIN  structure,  with  an  exponential 
increase  and  decrease  in  voltage.  The  dc  values  of  current  and  voltage  for  this  calculation 
are  V(DC)  =2.0v,  J(DC)  =4.22x10*  AJ cm2. 
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Figure  6.3b  Time  sequence  of  potential  proiles  corresponding  to  the  time  dependent 
calculation  of  figure  6.3a. 
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Figure  6.3c  Time  sequence  of  density  proiles  corresponding  to  the  time  dependent 
calculation  of  figure  6.3a. 
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Figure  6.3d  Time  sequence  of  gamma  valley  velocity  profiles  corresponding  to  the  time 
dependent  calculation  of  figure  6.3a 
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Figure  6.3e  Time  sequence  of  gamma  valley  temperature  profiles  corresponding  to  the  time 
dependent  calculation  of  figure  6.3  a. 
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The  comparative  calculations  of  figures  6.2  and  6.3  demonstrated  that  pulses  of  the  order 
of  five  picoseconds  in  duration  were  not  sufficient  to  elucidate  information  concerning 
velocity  overshoot,  but  that  intrapulse  structure  of  the  order  of  two  picoseconds  could 
emphasize  overshoot  contributions.  An  interesting  feature  of  the  calculation  was  that 
there  was  a  hint  that  the  structure  of  the  current  oscillation  could  provide  detailed 
differences  hbout  the  rate  of  intervalley  transfer  for  carriers  transferring  from  the  gamma 
to  L  valleys  and/or  from  the  L  to  the  gamma  valleys.  Since  it  is  known  that  the  most 
significant  effects  of  this  transfer  manifest  themselves  when  transfer  is  initiated,  the 
calculation  of  figure  63  was  repeated  but  with  a  lower  DC  bias  voltage  of  one  volt. 

The  calculations  are  shown  in  figure  6.4,  which  should  be  compared  to  figure  63.  In 
comparing  we  note  several  points.  There  are:  (i)  two  displacement  current  contributions, 
(ii)  higher  steady  state  values  of  current,  and  (iii)  the  differences  between  the  steady  state 
values  on  the  up  and  downswings  of  voltage.  But  unlike  the  calculations  of  figure  6.3,  the 
steady  state  current  differences  are  more  apparent  in  this  calculation.  The  details  of  this 
calculation  follow. 

Figure  6.4b  displays  the  distribution  of  voltage  across  the  structure  as  a  function  of  time. 
There  is  still  a  small  retarding  field  at  the  upstream  N +  region,  and  note  the  strong 
variation  in  potential  at  the  downstream  boundary.  The  profiles  of  electron  density, 
displayed  in  figure  6.4c  show  electron  transfer  near  the  anode  N  +  region,  but  significantly 
less  than  the  figure  63c  calculation.  One  the  basis  of  these  profiles,  we  note  that  at  the 
end  of  three  picoseconds  there  was  a  gradual  transfer  of  carriers  from  the  gamma  to  L 
valley,  while  in  the  fourth  picosecond  there  was  a  rapid  return  of  carriers  from  the  L  to  the 
gamma  valley.  The  gamma  valley  velocity  profiles  show  a  more  uniform  distribution  of 
velocity  across  the  entire  undoped  region,  with  values  near  4xl07  cm/sec,  with  higher 
values  near  the  anode  occurring  later  in  the  cycle  when  the  transfer  is  greatest  to  the  L 
valley.  Because  of  the  slow  electron  transfer  as  the  voltage  increases,  the  corresponding 
high  velocities  on  the  upswing;  and  the  rapid  transfer  of  carriers  from  the  L  to  gamma 
valleys  on  the  downswing,  the  steady  current  on  the  upswing  is  considerably  larger  than 
the  current  on  the  downswing.  Note:  the  electron  temperature,  figure  6.4e,  shows  the 
skewed  values  near  the  anode  with  most  of  the  values  hovering  around  the  point  where 
electron  transfer  should  occur.  The  comparison  of  figures  6.3  and  6.4  suggest  that  if  the 
relaxed  current  differences  during  one  cycle  can  be  measured  we  will  have  an  unequivocal 
measurement  of  nonequilibrium  transients  in  semiconductors. 

There  are  several  other  points  that  must  be  addressed,  one  concerned  with  the 
dependence  of  capacitive  or  displacement  current  effects  on  doping,  and  the  second  being 
the  effect  of  dealing  with  a  voltage  source  whose  period  is  significantly  shorter  than  the 
energy  relaxation  time.  Consider  first  displacement  current  effects. 
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Figure  6.4a  Voltage  and  current  profiles  for  the  nominal  N1N  structure,  with  an  exponential 
increase  and  decrease  in  voltage.  The  dc  values  of  current  and  voltage  for  this  calculation 
are  V(DC)=1.0v,  J(DC)  =3.57x10*  Alcm2. _ _ 
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Figure  6.4b  Time  sequence  of  potential  proiles  corresponding  to  the  time  dependent 
calculation  of  figure  6.4a. 
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Figure  6.4c  Time  sequence  of  density  proiles  corresponding  to  the  time  dependent  calculation 
of  figure  6.4a. 
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Figure  6.4d  Time  sequence  of  gamma  valley  velocity  profiles  corresponding  to  the  time 
dependent  calculation  of  figure  6.4a, 
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The  calculations  of  figure  6.4  were  repeated  but  with  the  doping  at  the  N  +  regions 
increased  to  101#/cms.  It  was  anticipated  that  the  displacement  current  oscillations, 
which  reflect  the  dielectric  relaxation  time  would  end  during  a  shorter  time  duration.  This 
occurred  and  is  shown  in  figure  6.5.  The  charge  distributions,  potential,  velocity  and 
temperature  are  strikingly  similar  to  the  distributions  for  the  lower  10 1 7 /cm3  calculation 
and  are  not  repeated. 

The  situation  when  the  figure  63  calculations  are  repeated  but  for  a  period  of  one 
picosecond  was  considered  next,  with  the  time  dependent  results  shown  in  figure  6.5.  At 
the  end  of  8  picoseconds  the  structure  did  not  reach  steady  state  equilibrium.  Indeed  the 
large  current  swings  are  essentially  displacement  current  effects,  and  while  velocity 
overshoot  is  occurring  at  all  times  during  this  oscillation,  it  would  be  difficult  to 
distinguish  this  contribution  from  the  displacement  current  contributions.  Indeed  it 
would  appear  that  the  measurement  of  velocity  overshoot  should  concentrate  on 
determining  the  changes  that  occur  in  going  from  steady  state  equilibrium  to  the  onset  on 
nonequilibrium. 

The  broad  conclusion  that  can  be  drawn  from  these  studies  is  that  nonequilibrium  transient 
effects  can  be  measured;  the  measurement  will  not  yield  overshoot  but  will  provide  direct 
confirmation  of  the  time  constants  associated  with  intervalley  transfer,  which  is  at  the  core  of 
nonequilibrium  phenomena  in  III-V materials. 
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Figure  6.5  Voltage  and  current  profiles  for  the  NIN  structure  with  more  heavily  doped  end 
regions,  with  an  exponential  increase  and  decrease  in  voltage.  The  dc  values  of  current  and 
voltage  for  this  calculation  are  V(DC)  =  l.Ov,  J(DC)  =3.63x10*  A/cm2. 
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Figure  6.6  Voltage  and  current  profiles  for  the  nominal  NIN  structure,  with  an  exponential 
increase  and  decrease  in  voltage.  The  period  of  this  oscllation  is  one  picosecond.  The  dc 
values  of  current  and  voltage  for  this  calculation  are  V(DC)=1.0v, 
J(DC)  =3.57x10*  A/cm2. 
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7.  Summary 

This  study  focused  on  the  implementation  of  numerical  algorithms  for  examining  the 
transient  behavior  of  electronic  devices  and  also  suggesting  experiments  to  elucidate 
velocity  overshoot. 

During  the  course  of  the  study  we  also  implemented  a  Monte  Carlo  algorithm  and 
enhanced  the  existing  moment  algorithm  to  treat  transients  more  accurately,  and  to 
include  two  dimensional  scattering  events. 

The  highlights  of  the  study  included: 

(1)  developing  new  scattering  integrals  for  two  dimensional  transport; 

(2)  establishing  the  presence  of  terahertz  frequency  small  signal  charge  density  waves 
forming  in  the  channel  of  the  two  dimensional  electron  gas;  and 

(3)  determining  that  measurements  could  be  made  that  would  expose  the  relative  time 
scale  differences  between  energy  and  momentum  relaxation  effects-thereby  providing 
evidence  for  nonequilibrium  velocity  overshoot. 
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